7. ADI+IFS: PSF subtraction and forward modeling

Author: Valentin Christiaens
Suitable for VIP v1.2.2 onwards (Warning this is a different version requirement than other tutorials)
Last update: 2022/04/13

Table of contents

The purpose of this tutorial is to show: - how to inject a planet spectrum in a 4D (IFS+ADI) datacube; - how to post-process 4D (IFS+ADI) datacubes; - how to retrieve the parameters of a planet present in a 4D (IFS+ADI) datacube (astro- and spectrometry).

Note that the model spectra (for the star and planet) used in this notebook are provided after resampling for the purpose of highlighting the capabilities of VIP only (i.e. to avoid having to install additional requirements to be able to run this tutorial). If you are interested in a more complete planet spectrum injection procedure, check out this Exoplanet Data Challenge Phase 2 tutorial, which also requires the following external packages: `special <https://github.com/VChristiaens/special>`__ (for model resampling) and `eidc2 <https://github.com/exoplanet-imaging-challenge/phase2>`__ (for better flux normalization before spectrum injection).

Let’s first load all packages needed in this tutorial:

[1]:

import astropy.constants as c
import glob
from hciplot import plot_frames, plot_cubes
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd

#from special.model_resampling import resample_model
#from special.utils_spec import find_nearest

from vip_hci.config import VLT_SPHERE_IFS
from vip_hci.fits import open_fits, write_fits
from vip_hci.fm import (cube_inject_companions, firstguess, mcmc_negfc_sampling, normalize_psf,
                        show_walk_plot, show_corner_plot, confidence)
from vip_hci.metrics import snr, significance
from vip_hci.psfsub import median_sub, pca, pca_annular, pca_annulus, pca_grid
from vip_hci.preproc import find_scal_vector, frame_rescaling
from vip_hci.var import frame_center, dist, mask_circle

7.1. Loading and visualizing the data

In the ‘datasets’ folder of the VIP_extras repository you can find a toy SPHERE/IFS coronagraphic cube acquired in pupil-stabilized mode on the source HIP39826 (a star with no reported directly imaged companion). The folder also contains the associated non-coronagraphic point spread function (PSF), wavelength vector, parallactic angles, and airmass of the source.

Let’s now load the data. Note that more info on opening and visualizing fits files with VIP in general is available in the first VIP tutorial.

[2]:

path = '../datasets/'

cubename = path+'sphere_ifs_HIP39826_cube.fits'
psfname = path+'sphere_ifs_HIP39826_psf.fits'
lbdaname = path+'sphere_ifs_HIP39826_wl.fits'
angname = path+'sphere_ifs_HIP39826_pa.fits'
Xname = path+'sphere_ifs_HIP39826_airmass.fits'

cube = open_fits(cubename)
psf = open_fits(psfname)
lbda = open_fits(lbdaname)
derot_angles = -1*open_fits(angname) # note: SPHERE DRH returns parallactic angles, while VIP routines take derotation angles as input
X_sci = open_fits(Xname)

nch, nz, ny, nx = cube.shape

Fits HDU-0 data successfully loaded. Data shape: (39, 55, 107, 107)
Fits HDU-0 data successfully loaded. Data shape: (39, 29, 29)
Fits HDU-0 data successfully loaded. Data shape: (39,)
Fits HDU-0 data successfully loaded. Data shape: (55,)
Fits HDU-0 data successfully loaded. Data shape: (55,)

Each IFS spectral cube consists of 39 monochromatic images spread in wavelengths between the Y and J bands (‘YJ’ mode) or Y and H bands (‘YJH’ mode). Here the IFS+ADI cube contains 55 such spectral cubes combined into a single master cube.

The master spectral cube is already centered, trimmed from bad frames, and mostly corrected from bad pixels (check Tutorial 2 for example uses).

Let’s inspect the first wavelength using hciplot.plot_cubes (feel free to set the backend to ‘bokeh’ to read pixel values interactively):

[3]:

plot_cubes(cube)#, backend='bokeh')

:Dataset   [x,y,time,lambda]   (flux)
:Cube_shape     [107, 107, 55, 39]

[3]:

We notice that the flux varies and the PSF size increases, both as a function of wavelength.

Now let’s inspect the PSF cube (one frame per wavelength), which will be used for the actual injection of a fake planet:

[4]:

plot_cubes(psf)

:Dataset   [x,y,time]   (flux)
:Cube_shape     [29, 29, 39]

[4]:

As we can see, the flux varies as a function of wavelength. This reflects the intrinsic spectrum of the star but also the instrumental transmission and atmospheric absorption bands (e.g. associated to the water molecule) affecting some wavelengths more than others. The flux will be normalized in the next section, before injection.

Now let’s check the parallactic angle and airmass variations:

[5]:

%matplotlib inline
fig, axes = plt.subplots(1,2, figsize = (14,5))
ax1, ax2 = axes
ax1.plot(range(1,nz+1), derot_angles)
ax1.set_xlabel("Cube index")
ax1.set_ylabel("Parallactic angle (deg)")
ax2.plot(range(1,nz+1), X_sci)
ax2.set_xlabel("Cube index")
ax2.set_ylabel("Airmass")
plt.show()
../_images/tutorials_07_ifs_psfsub_fm_planets_18_0.png

Note the sudden variation in parallactic angle at the end of the sequence is because we trimmed bad frames from the sequence (all were located towards the end of the observation).

Let’s define the pixel scale for IFS - note that this is optional since the planet coordinates will be given in pixels. It will still be useful to have meaningful radial separations expressed in arcsec.

[6]:

plsc = VLT_SPHERE_IFS['plsc']
print(plsc, "arcsec/px")

0.0074 arcsec/px

Let’s also define the mean spectral resolution of the instrument (which will be used for convolution of model spectra):

[7]:

spec_res = 50.

Let’s now load and visualize the coronagraphic radial transmission profile:

[8]:

trans = open_fits(path+"../datasets/sphere_ifs_YJ_APLC_CoroTransmission.fits")

Fits HDU-0 data successfully loaded. Data shape: (40, 23)

[9]:

last_r = 20
plt.figure(figsize=(10,6))
plt.plot(trans[0,:last_r]*1000, trans[1,:last_r], 'b', label='ch 1')
plt.plot(trans[0,:last_r]*1000, trans[10,:last_r], 'c', label='ch 10')
plt.plot(trans[0,:last_r]*1000, trans[20,:last_r], 'y', label='ch 20')
plt.plot(trans[0,:last_r]*1000, trans[30,:last_r], 'm', label='ch 30')
plt.plot(trans[0,:last_r]*1000, trans[-1,:last_r], 'r', label='ch 39')
plt.xlabel("Radius (mas)")
plt.ylabel("Transmission")
plt.legend()

[9]:

<matplotlib.legend.Legend at 0x7f8d3020afd0>
../_images/tutorials_07_ifs_psfsub_fm_planets_26_1.png

EDIT: this is not the most up-to-date SPHERE coronagraphic transmission curve. Please contact relevant SPHERE colleagues or check the ESO website for a newer version. Do not use the above curve for your research.

7.2. Planet parameters and spectra

Now let’s set the ground truth coordinates and average contrast (over all wavelengths) at which two planets will be injected in this cube:

[10]:

xy_b = (58.6,37.7)
contrast_b = 1e-3

xy_c = (83.8,29.9)
contrast_c = 2e-4

Let’s convert cartesian to polar coordinates:

[11]:

cy, cx = frame_center(cube)
r_b = dist(cy, cx, xy_b[1], xy_b[0])
r_c = dist(cy, cx, xy_c[1], xy_c[0])

theta_b = np.rad2deg(np.arctan2(xy_b[1]-cy, xy_b[0]-cx))
theta_c = np.rad2deg(np.arctan2(xy_c[1]-cy, xy_c[0]-cx))

[12]:

print(cy, cx)
print(r_b, theta_b)
print(r_c, theta_c)

53 53
16.292636373527763 -69.89671372878468
38.5 -36.86989764584403

Now let’s consider two fiducial model spectra for the injection of the 2 planets: a flat spectrum and a 1000K blackbody model (note that this exercise is purely academic).

[13]:

def blackbody(lbda, T):
    fac = 2*c.h.value*(c.c.value**2)/(np.power(lbda*1e-6,5))
    div = (1/(np.exp((c.h.value*c.c.value)/((lbda*1e-6)*c.k_B.value*T))-1))
    # convert from W m-3 Sr-1 to W m-2 mu-1 Sr-1
    conv = 1e-6
    return fac*div*conv

[14]:

fluxes_b = np.ones_like(lbda)
fluxes_c = blackbody(lbda, 1000)

spec_b = np.array([lbda, fluxes_b])
spec_c = np.array([lbda, fluxes_c])

[15]:

plt.plot(lbda, fluxes_b/np.amax(fluxes_b), 'k', label='flat')
plt.plot(lbda, fluxes_c/np.amax(fluxes_c), 'r', label='1000 K BB')
plt.xlabel(r"$\lambda (\mu$m)")
plt.ylabel(r"$ F_{\lambda}$ (arbitrary units)")
plt.legend()

[15]:

<matplotlib.legend.Legend at 0x7f8d425d76a0>
../_images/tutorials_07_ifs_psfsub_fm_planets_37_1.png

7.3. Flux normalization

7.3.1. Non-coronagraphic PSF normalization

Let’s measure the flux and FWHM of the non-coronagraphic PSF and normalize its flux to 1 in a 1-FWHM radius aperture (all at once), in all channels, using VIP’s fm/normalize_psf function. We use the debug option to check that each 2D Gaussian fit worked properly:

[16]:

%matplotlib inline
psfn, flux_st, fwhm = normalize_psf(psf, fwhm='fit', full_output=True, debug=True)
../_images/tutorials_07_ifs_psfsub_fm_planets_41_0.png 
FWHM_y = 6.249918590625379
FWHM_x = 5.955174570821388

centroid y = 13.912712804514848
centroid x = 13.89930852239608
centroid y subim = 13.912712804514848
centroid x subim = 13.89930852239608

amplitude = 636.3209207145087
theta = -40.20327494950595
../_images/tutorials_07_ifs_psfsub_fm_planets_41_2.png 
FWHM_y = 6.115031733427154
FWHM_x = 5.83909291219867

centroid y = 14.027611189747015
centroid x = 13.902123414264723
centroid y subim = 14.027611189747015
centroid x subim = 13.902123414264723

amplitude = 1054.2949916064315
theta = -37.536071255420595
../_images/tutorials_07_ifs_psfsub_fm_planets_41_4.png 
FWHM_y = 6.017536940978627
FWHM_x = 5.754055056927457

centroid y = 14.064675014828953
centroid x = 13.92189940782825
centroid y subim = 14.064675014828953
centroid x subim = 13.92189940782825

amplitude = 1404.7946951267922
theta = -32.753406960100804
../_images/tutorials_07_ifs_psfsub_fm_planets_41_6.png 
FWHM_y = 5.9918781552060665
FWHM_x = 5.682647009422151

centroid y = 14.065184394206616
centroid x = 13.946038597175635
centroid y subim = 14.065184394206616
centroid x subim = 13.946038597175635

amplitude = 1705.2716467832417
theta = -37.38382831177306
../_images/tutorials_07_ifs_psfsub_fm_planets_41_8.png 
FWHM_y = 6.025210456837181
FWHM_x = 5.719302651407007

centroid y = 14.075784813017359
centroid x = 13.955473078460217
centroid y subim = 14.075784813017359
centroid x subim = 13.955473078460217

amplitude = 1975.512023782089
theta = -39.227168101522786
../_images/tutorials_07_ifs_psfsub_fm_planets_41_10.png 
FWHM_y = 5.985037638840228
FWHM_x = 5.71990434228688

centroid y = 14.090563585852108
centroid x = 13.954901979145388
centroid y subim = 14.090563585852108
centroid x subim = 13.954901979145388

amplitude = 2143.996699880163
theta = -37.89243906945627
../_images/tutorials_07_ifs_psfsub_fm_planets_41_12.png 
FWHM_y = 6.037449649434506
FWHM_x = 5.753733367061647

centroid y = 14.109510955089752
centroid x = 13.980317492176528
centroid y subim = 14.109510955089752
centroid x subim = 13.980317492176528

amplitude = 2165.8668020756913
theta = -38.98816537806606
../_images/tutorials_07_ifs_psfsub_fm_planets_41_14.png 
FWHM_y = 6.010425979658148
FWHM_x = 5.749988204487918

centroid y = 14.091399025168501
centroid x = 13.986593684665712
centroid y subim = 14.091399025168501
centroid x subim = 13.986593684665712

amplitude = 2237.4526741829254
theta = -42.08255560945647
../_images/tutorials_07_ifs_psfsub_fm_planets_41_16.png 
FWHM_y = 6.005151597186863
FWHM_x = 5.766189137359492

centroid y = 14.090254937281157
centroid x = 14.000261828564378
centroid y subim = 14.090254937281157
centroid x subim = 14.000261828564378

amplitude = 2355.1073353636825
theta = -38.12872688398152
../_images/tutorials_07_ifs_psfsub_fm_planets_41_18.png 
FWHM_y = 5.968097471112488
FWHM_x = 5.732989798766061

centroid y = 14.102026780112409
centroid x = 14.00192256597961
centroid y subim = 14.102026780112409
centroid x subim = 14.00192256597961

amplitude = 2523.955221661578
theta = -36.16671879464335
../_images/tutorials_07_ifs_psfsub_fm_planets_41_20.png 
FWHM_y = 6.201045253062894
FWHM_x = 5.952529125946643

centroid y = 14.09561886484433
centroid x = 14.011238141279097
centroid y subim = 14.09561886484433
centroid x subim = 14.011238141279097

amplitude = 2483.278573207968
theta = -40.01453306537682
../_images/tutorials_07_ifs_psfsub_fm_planets_41_22.png 
FWHM_y = 6.1167439995646244
FWHM_x = 5.839271620946151

centroid y = 14.100847253490924
centroid x = 14.00962801008153
centroid y subim = 14.100847253490924
centroid x subim = 14.00962801008153

amplitude = 2649.287384069412
theta = -42.97211525180409
../_images/tutorials_07_ifs_psfsub_fm_planets_41_24.png 
FWHM_y = 6.037056836628313
FWHM_x = 5.786392745942098

centroid y = 14.110142788418317
centroid x = 14.004784235107909
centroid y subim = 14.110142788418317
centroid x subim = 14.004784235107909

amplitude = 2770.0240286010007
theta = -39.87700223659835
../_images/tutorials_07_ifs_psfsub_fm_planets_41_26.png 
FWHM_y = 6.066318705156946
FWHM_x = 5.822724547170392

centroid y = 14.099317823625379
centroid x = 14.003150239110157
centroid y subim = 14.099317823625379
centroid x subim = 14.003150239110157

amplitude = 2830.2548554359573
theta = -47.81201829549842
../_images/tutorials_07_ifs_psfsub_fm_planets_41_28.png 
FWHM_y = 6.076253484626576
FWHM_x = 5.8064868728535695

centroid y = 14.112521461380666
centroid x = 14.034377198493269
centroid y subim = 14.112521461380666
centroid x subim = 14.034377198493269

amplitude = 2803.405353990872
theta = -1130.8992713589835
../_images/tutorials_07_ifs_psfsub_fm_planets_41_30.png 
FWHM_y = 5.78340826410417
FWHM_x = 6.094768878919073

centroid y = 14.151451481603505
centroid x = 14.04623588538711
centroid y subim = 14.151451481603505
centroid x subim = 14.04623588538711

amplitude = 2553.3202960394237
theta = -1046.1876651647306
../_images/tutorials_07_ifs_psfsub_fm_planets_41_32.png 
FWHM_y = 6.181608395179831
FWHM_x = 5.8193909503216625

centroid y = 14.206064453326714
centroid x = 14.027564826853363
centroid y subim = 14.206064453326714
centroid x subim = 14.027564826853363

amplitude = 2063.1990503463253
theta = -591.8971102346304
../_images/tutorials_07_ifs_psfsub_fm_planets_41_34.png 
FWHM_y = 5.751799173200897
FWHM_x = 6.1765747725042495

centroid y = 14.199498957513265
centroid x = 14.014692732639789
centroid y subim = 14.199498957513265
centroid x subim = 14.014692732639789

amplitude = 1668.5839475163139
theta = 214.19574599466478
../_images/tutorials_07_ifs_psfsub_fm_planets_41_36.png 
FWHM_y = 5.797247572765799
FWHM_x = 6.140720682799337

centroid y = 14.09188009913005
centroid x = 13.97827848457926
centroid y subim = 14.09188009913005
centroid x subim = 13.97827848457926

amplitude = 1587.9799969490432
theta = 208.9336267041074
../_images/tutorials_07_ifs_psfsub_fm_planets_41_38.png 
FWHM_y = 6.1732867167382
FWHM_x = 5.8497702682633586

centroid y = 14.020275158168136
centroid x = 14.00256734343732
centroid y subim = 14.020275158168136
centroid x subim = 14.00256734343732

amplitude = 1805.5945893071935
theta = 117.55114923044374
../_images/tutorials_07_ifs_psfsub_fm_planets_41_40.png 
FWHM_y = 5.879233953058072
FWHM_x = 6.127221747592718

centroid y = 13.989753883923711
centroid x = 14.002807252670436
centroid y subim = 13.989753883923711
centroid x subim = 14.002807252670436

amplitude = 2198.1572616046824
theta = -320.21556686945866
../_images/tutorials_07_ifs_psfsub_fm_planets_41_42.png 
FWHM_y = 5.887171925622423
FWHM_x = 6.160164902722713

centroid y = 14.006124941703272
centroid x = 14.009936779665049
centroid y subim = 14.006124941703272
centroid x subim = 14.009936779665049

amplitude = 2628.4500069112014
theta = -321.3590169947588
../_images/tutorials_07_ifs_psfsub_fm_planets_41_44.png 
FWHM_y = 6.0983193806752976
FWHM_x = 5.823811537198284

centroid y = 14.04578078126952
centroid x = 14.01477101186026
centroid y subim = 14.04578078126952
centroid x subim = 14.01477101186026

amplitude = 2992.2717104226613
theta = -228.95610743164002
../_images/tutorials_07_ifs_psfsub_fm_planets_41_46.png 
FWHM_y = 6.088143522859679
FWHM_x = 5.8003209490757035

centroid y = 14.051247831826354
centroid x = 13.998093842074548
centroid y subim = 14.051247831826354
centroid x subim = 13.998093842074548

amplitude = 3191.2868184142203
theta = -409.93971951101645
../_images/tutorials_07_ifs_psfsub_fm_planets_41_48.png 
FWHM_y = 6.131166433296467
FWHM_x = 5.872212732534487

centroid y = 14.035337815431665
centroid x = 13.998171477359165
centroid y subim = 14.035337815431665
centroid x subim = 13.998171477359165

amplitude = 3243.165884766142
theta = -50.09403204938458
../_images/tutorials_07_ifs_psfsub_fm_planets_41_50.png 
FWHM_y = 6.1749426650386505
FWHM_x = 5.9932124659024835

centroid y = 14.039753518526705
centroid x = 13.991921001229356
centroid y subim = 14.039753518526705
centroid x subim = 13.991921001229356

amplitude = 3307.8703264494
theta = -60.778524959138785
../_images/tutorials_07_ifs_psfsub_fm_planets_41_52.png 
FWHM_y = 6.212946735034473
FWHM_x = 6.032862556063642

centroid y = 14.009742406995104
centroid x = 13.967991908973474
centroid y subim = 14.009742406995104
centroid x subim = 13.967991908973474

amplitude = 3436.0292600146126
theta = -233.074343782064
../_images/tutorials_07_ifs_psfsub_fm_planets_41_54.png 
FWHM_y = 6.202834709617493
FWHM_x = 5.991454580377878

centroid y = 14.020014512640392
centroid x = 13.952064722414017
centroid y subim = 14.020014512640392
centroid x subim = 13.952064722414017

amplitude = 3546.656042036664
theta = -45.08873953984778
../_images/tutorials_07_ifs_psfsub_fm_planets_41_56.png 
FWHM_y = 6.188582125733534
FWHM_x = 5.985402654635661

centroid y = 14.012885928141898
centroid x = 13.933314914222644
centroid y subim = 14.012885928141898
centroid x subim = 13.933314914222644

amplitude = 3607.1363107553057
theta = -45.542552785379904
../_images/tutorials_07_ifs_psfsub_fm_planets_41_58.png 
FWHM_y = 6.244755229551024
FWHM_x = 6.022172272876057

centroid y = 14.00590897381181
centroid x = 13.928682611420799
centroid y subim = 14.00590897381181
centroid x subim = 13.928682611420799

amplitude = 3616.1909624995233
theta = -32.925763043167755
../_images/tutorials_07_ifs_psfsub_fm_planets_41_60.png 
FWHM_y = 6.27304172170033
FWHM_x = 6.080609549792838

centroid y = 13.982299985250767
centroid x = 13.934886096629457
centroid y subim = 13.982299985250767
centroid x subim = 13.934886096629457

amplitude = 3579.548593229253
theta = -38.810946878862445
../_images/tutorials_07_ifs_psfsub_fm_planets_41_62.png 
FWHM_y = 6.272897562677739
FWHM_x = 6.150604077395562

centroid y = 13.949986869106837
centroid x = 13.940065010214068
centroid y subim = 13.949986869106837
centroid x subim = 13.940065010214068

amplitude = 3535.36267540835
theta = -1487.052254957365
../_images/tutorials_07_ifs_psfsub_fm_planets_41_64.png 
FWHM_y = 6.169653547163166
FWHM_x = 6.268150600253964

centroid y = 13.931026676902302
centroid x = 13.920174694388791
centroid y subim = 13.931026676902302
centroid x subim = 13.920174694388791

amplitude = 3562.252494864886
theta = 32.842817255455444
../_images/tutorials_07_ifs_psfsub_fm_planets_41_66.png 
FWHM_y = 6.274537554500622
FWHM_x = 6.157837787311417

centroid y = 13.934250084708902
centroid x = 13.92169790717976
centroid y subim = 13.934250084708902
centroid x subim = 13.92169790717976

amplitude = 3655.664070910777
theta = 129.0032036755185
../_images/tutorials_07_ifs_psfsub_fm_planets_41_68.png 
FWHM_y = 6.311436197874585
FWHM_x = 6.135315847492893

centroid y = 13.945184682040301
centroid x = 13.90455843917349
centroid y subim = 13.945184682040301
centroid x subim = 13.90455843917349

amplitude = 3713.222558796391
theta = -41.34373753210502
../_images/tutorials_07_ifs_psfsub_fm_planets_41_70.png 
FWHM_y = 6.330588326155178
FWHM_x = 6.183593806815872

centroid y = 13.97269314976903
centroid x = 13.886442080617911
centroid y subim = 13.97269314976903
centroid x subim = 13.886442080617911

amplitude = 3603.9831790818416
theta = -39.62090236411358
../_images/tutorials_07_ifs_psfsub_fm_planets_41_72.png 
FWHM_y = 6.184001484216053
FWHM_x = 6.343494966712782

centroid y = 13.949799220444797
centroid x = 13.882979812466415
centroid y subim = 13.949799220444797
centroid x subim = 13.882979812466415

amplitude = 3362.492517413177
theta = -131.56408968416724
../_images/tutorials_07_ifs_psfsub_fm_planets_41_74.png 
FWHM_y = 6.104938340342794
FWHM_x = 6.326105574311393

centroid y = 13.963528380019108
centroid x = 13.869285619906863
centroid y subim = 13.963528380019108
centroid x subim = 13.869285619906863

amplitude = 3062.496165602359
theta = -498.25201487357776
../_images/tutorials_07_ifs_psfsub_fm_planets_41_76.png 
FWHM_y = 6.281872522029886
FWHM_x = 6.106896271562541

centroid y = 14.035165570902214
centroid x = 13.854666632379166
centroid y subim = 14.035165570902214
centroid x subim = 13.854666632379166

amplitude = 2583.840814837761
theta = -35.5994403427372
Mean FWHM per channel:
[6.103 5.977 5.886 5.837 5.872 5.852 5.896 5.88  5.886 5.851 6.077 5.978
 5.912 5.945 5.941 5.939 6.    5.964 5.969 6.012 6.003 6.024 5.961 5.944
 6.002 6.084 6.123 6.097 6.087 6.133 6.177 6.212 6.219 6.216 6.223 6.257
 6.264 6.216 6.194]
Flux in 1xFWHM aperture:
[13397.649 21373.517 27661.433 33057.763 38754.981 41781.661 42805.645
 44021.68  46395.703 49082.614 52371.307 53918.696 55093.806 56930.742
 56350.541 51298.303 42189.297 33733.489 32181.659 37142.073 45223.932
 54394.967 60556.863 64201.073 66477.378 69675.456 73203.718 74978.343
 76002.683 77311.678 77583.749 77448.17  78182.432 80187.787 81616.948
 80114.272 74899.504 67091.76  56185.008]

Now let’s plot the measured FWHM and stellar flux:

[17]:

fig, axes = plt.subplots(1,2, figsize = (14,5))
ax1, ax2 = axes
ax1.plot(lbda, fwhm, 'bo')
ax1.set_xlabel(r"Wavelength ($\mu$m)")
ax1.set_ylabel("FWHM")
ax2.plot(lbda, flux_st, 'ro')
ax2.set_xlabel(r"Wavelength ($\mu$m)")
ax2.set_ylabel("Stellar flux (ADUs)")
plt.show()
../_images/tutorials_07_ifs_psfsub_fm_planets_43_0.png

Fig. 3.1: Measured FWHM and stellar flux in the non-coronagraphic PSF cube.

[18]:

plot_cubes(psfn, grid=True)

:Dataset   [x,y,time]   (flux)
:Cube_shape     [29, 29, 39]

[18]:

Now let’s save the measured spectrum of the star in the format we will need for the injection function:

[19]:

spec_star = np.array([lbda, flux_st])

7.3.2. Stellar SED model

The spectral type of the star (HIP39826) is K7V. For the purpose of this tutorial, we simply load a BT-NextGen SED model that has already be resampled considering the spectral resolution and sampling of SPHERE-IFS:

[20]:

spec_star_mod = open_fits("../datasets/sphere_ifs_HIP39826_SED_model.fits")

plt.plot(spec_star_mod[0], spec_star_mod[0]*spec_star_mod[1], 'k', label='Model SED for the star')
plt.ylabel(r"$\lambda F_{\lambda}$ (arbitrary units)")
plt.legend()

Fits HDU-0 data successfully loaded. Data shape: (2, 39)

[20]:

<matplotlib.legend.Legend at 0x7f8d52a15460>
../_images/tutorials_07_ifs_psfsub_fm_planets_50_2.png

7.3.3. Airmass

Airmass variation can also be taken into consideration during the injection of the planet spectra (more specifically we consider the ratio between individual airmass values and the median airmass over the observation):

[21]:

med_X = np.median(X_sci)
X_fac = np.exp(-(X_sci-med_X))

Let’s plot the factors to be applied to correct for airmass (w.r.t median airmass):

[22]:

%matplotlib inline
plt.plot(range(1,nz+1), X_fac, 'k')
plt.xlabel("Cube index")
plt.ylabel("X factor")

[22]:

Text(0, 0.5, 'X factor')
../_images/tutorials_07_ifs_psfsub_fm_planets_55_1.png

7.3.4. Final fluxes

Let’s now calculate the fluxes at which the normalized PSF will actually be injected in each frame of each spectral cube.

Considering the effect of instrumental and atmospheric transmission, we have:

[23]:

flux_b_trans = fluxes_b*(flux_st/spec_star_mod[1])
flux_c_trans = fluxes_c*(flux_st/spec_star_mod[1])

Considering the actual mean contrast we want the planet to be injected, with respect to the star, we have:

[24]:

flux_b_scal = flux_b_trans*contrast_b*np.mean(flux_st)/np.mean(flux_b_trans)
flux_c_scal = flux_c_trans*contrast_c*np.mean(flux_st)/np.mean(flux_c_trans)

Finally let’s consider the airmass:

[25]:

final_fluxes_b = np.outer(flux_b_scal, X_fac)
final_fluxes_c = np.outer(flux_c_scal, X_fac)

Note that the effect of coronagraphic transmission is directly taken into account in the injection function (next section).

7.4. Injection

The actual injection can be simply performed with the cube_inject_companions function, which can handle 4D input cubes. Let’s inject the first planet first:

[26]:

cube_fc = cube_inject_companions(cube, psfn, derot_angles, flevel=final_fluxes_b,
                                 plsc=plsc, rad_dists=[r_b], n_branches=1,
                                 theta=theta_b, transmission=trans, verbose=True)

*** Processing spectral channel 1/39 ***
*** Processing spectral channel 2/39 ***
*** Processing spectral channel 3/39 ***
*** Processing spectral channel 4/39 ***
*** Processing spectral channel 5/39 ***
*** Processing spectral channel 6/39 ***
*** Processing spectral channel 7/39 ***
*** Processing spectral channel 8/39 ***
*** Processing spectral channel 9/39 ***
*** Processing spectral channel 10/39 ***
*** Processing spectral channel 11/39 ***
*** Processing spectral channel 12/39 ***
*** Processing spectral channel 13/39 ***
*** Processing spectral channel 14/39 ***
*** Processing spectral channel 15/39 ***
*** Processing spectral channel 16/39 ***
*** Processing spectral channel 17/39 ***
*** Processing spectral channel 18/39 ***
*** Processing spectral channel 19/39 ***
*** Processing spectral channel 20/39 ***
*** Processing spectral channel 21/39 ***
*** Processing spectral channel 22/39 ***
*** Processing spectral channel 23/39 ***
*** Processing spectral channel 24/39 ***
*** Processing spectral channel 25/39 ***
*** Processing spectral channel 26/39 ***
*** Processing spectral channel 27/39 ***
*** Processing spectral channel 28/39 ***
*** Processing spectral channel 29/39 ***
*** Processing spectral channel 30/39 ***
*** Processing spectral channel 31/39 ***
*** Processing spectral channel 32/39 ***
*** Processing spectral channel 33/39 ***
*** Processing spectral channel 34/39 ***
*** Processing spectral channel 35/39 ***
*** Processing spectral channel 36/39 ***
*** Processing spectral channel 37/39 ***
*** Processing spectral channel 38/39 ***
*** Processing spectral channel 39/39 ***

Let’s now inject the second planet in the same cube:

[27]:

cube_fc = cube_inject_companions(cube_fc, psfn, derot_angles, flevel=final_fluxes_c,
                                 plsc=plsc, rad_dists=[r_c], n_branches=1,
                                 theta=theta_c, transmission=trans, verbose=False)

Let’s inspect the cube with the inject companions:

[28]:

plot_cubes(cube_fc)

:Dataset   [x,y,time,lambda]   (flux)
:Cube_shape     [107, 107, 55, 39]

[28]:

Go to the top

7.5. Post-processing of IFS data

Let’s now try to recover the 2 injected companions using different post-processing algorithms.

7.5.1. Optimal scaling factors

In order to leverage the radial expansion of the stellar PSF with wavelength through spectral differential imaging, let’s first compute the optimal scaling factors between the spectral channels. For a perfect AO correction and flat stellar spectrum, one would expect the scaling factors to be (exactly) inversely proportional to the wavelength. Reality deviates a bit from these assumptions but this still makes for a good first guess:

[29]:

th_scal = lbda[-1]/lbda
th_scal

[29]:

array([1.3877398, 1.375682 , 1.3635113, 1.3512615, 1.3389585, 1.3266245,
       1.3142908, 1.3019674, 1.2896852, 1.2774729, 1.2653329, 1.2532798,
       1.2413497, 1.2295421, 1.217879 , 1.2063586, 1.1950108, 1.1838219,
       1.1728301, 1.1620203, 1.151418 , 1.1410172, 1.1308315, 1.1208638,
       1.1111166, 1.1016016, 1.092311 , 1.0832641, 1.0744526, 1.0658859,
       1.0575558, 1.0494791, 1.0416555, 1.0340765, 1.0267497, 1.0196825,
       1.0128585, 1.0063006, 1.       ], dtype=float32)

Let’s use preproc.find_scal_vector to refine this estimate. This function will find the optimal scaling factors (both in terms of geometric scaling and flux scaling) between each spectral channel and the last one such that their subtraction yields minimal residuals. For coronagraphic observations, it is recommended to provide a binary mask where the residuals will be calculated which avoids the area under and at the edge of the coronagraphic mask for best performances. Here we use the mean of the 4D cube along the temporal direction to find the factors:

[30]:

mask = np.ones_like(cube[0,0])
mask = mask_circle(mask, int(0.16/plsc))

opt_scal_mean, opt_flux = find_scal_vector(np.mean(cube,axis=1), lbda, flux_st,
                                           mask=mask, nfp=2, fm="stddev")

Although a bit slower, another approach would be to calculate the scaling factors for each individual spectral cube, and use the median values:

[31]:

opt_scals = np.zeros([cube.shape[1],cube.shape[0]])
opt_fluxes = np.zeros([cube.shape[1],cube.shape[0]])
for i in range(cube.shape[1]):
    opt_scals[i], opt_fluxes[i] = find_scal_vector(cube[:,i], lbda, flux_st, mask=mask, nfp=2, fm="stddev")

[32]:

opt_scal_med = np.median(opt_scals,axis=0)
opt_flux_med = np.median(opt_fluxes,axis=0)

Let’s plot the theoretical and empirical scaling factors:

[33]:

plt.figure(figsize=(8,6))
for i in range(cube.shape[1]):
    if not i%10:
        label = 'cube #{}'.format(i)
    else:
        label = None
    plt.plot(lbda, opt_scals[i], alpha=0.3, label=label)
plt.plot(lbda, th_scal, 'k:', lw=2, label='Theoretical')
plt.plot(lbda, opt_scal_mean, 'k--', lw=2, label='Optimal (mean cube)')
plt.plot(lbda, np.median(opt_scals,axis=0), 'k-', lw=2, label='Optimal (median of indiv. curves)')
plt.xlabel(r"Wavelength ($\mu$m)")
plt.ylabel("Scaling factor")
plt.legend()

[33]:

<matplotlib.legend.Legend at 0x7f8d30460490>
../_images/tutorials_07_ifs_psfsub_fm_planets_84_1.png

Let’s visualize the actual residuals in the first spectral cube:

[34]:

res_cube_test_opt = np.zeros([cube.shape[0],cube.shape[-2],cube.shape[-1]])
for i in range(cube.shape[0]):
    res_cube_test_opt[i] = opt_fluxes[0,i]*frame_rescaling(cube[i,0], scale=opt_scal_med[i])-cube[-1,0]

[35]:

plot_cubes(res_cube_test_opt, vmin=-5)

:Dataset   [x,y,time]   (flux)
:Cube_shape     [107, 107, 39]

[35]:

Apart from the residual starlight at the edge of the coronagraphic mask (which does not expand radially), we see that the stellar halo is very well subtracted.

7.5.2. median-SDI

The most basic way to leverage the radial diversity present in IFS data is to perform median-SDI (i.e. the radial equivalent to median-ADI). For that just use the psfsub.median_sub function, with the sdi_only option on:

[36]:

mask_px = int(0.1/plsc)

med_SDI = median_sub(cube_fc, derot_angles, scale_list=opt_scal_med, sdi_only=True,
                     radius_int=mask_px, interp_zeros=True, mask_val=0)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:47:50
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
39 spectral channels per IFS frame
First median subtraction exploiting spectral variability
Running time:  0:00:17.046448
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
Median subtraction in the ADI fashion
Done derotating and combining
Running time:  0:00:19.216818
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

For better results, it is recommended to use both the geometrical and flux scaling factors found in the previous section:

[37]:

med_SDI_opt = median_sub(cube_fc, derot_angles, scale_list=opt_scal_med, flux_sc_list=opt_flux_med, sdi_only=True,
                         radius_int=mask_px, interp_zeros=True, mask_val=0)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:48:09
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
39 spectral channels per IFS frame
First median subtraction exploiting spectral variability
Running time:  0:00:17.223937
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
Median subtraction in the ADI fashion
Done derotating and combining
Running time:  0:00:19.351697
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s compare both reductions, showing circles at the ground truth locations where we injected the fake companions.

[38]:

%matplotlib inline
plot_frames((med_SDI, med_SDI_opt), circle=(xy_b,xy_c), vmin=-0.2, vmax=0.6)
../_images/tutorials_07_ifs_psfsub_fm_planets_95_0.png

In the first image, we notice the limitation of SDI at short angular separation. Since the stellar residuals at the edge of the coronagraphic mask do not expand radially as the rest of the stellar halo, we are left with a bright spurious ring, preventing the re-detection of fake planet b.

We see that the median-SDI reduction with scaled fluxes performs slightly better at subtracting stellar residuals near the edge of the coronagraphic mask, although fake planet b is still not recovered.

Fake planet c is redetected. Let’s compare its signal-to-noise ratio in both images:

[39]:

snr_c = snr(med_SDI, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
snr_c_sc = snr(med_SDI_opt, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c, snr_c_sc)

19.808257154238373 23.048229928839643

Overall, we can conclude that the scaled flux version of median-SDI performs better than the non-flux scaled version.

7.5.3. PCA-SDI

Let’s now run a PCA-SDI reduction. Note that this approach may be useful in datasets with azimuthally extended disc features (since ADI can filter such signals).

For PCA reductions, flux scaling factors are not necessary. This is for several reasons: i) the basis is orthonormal such that ; ii) an internal scaling is performed anyway when the principal components are projected onto each images; iii) flux scaling will typically increase the noise contribution before calculation of the principal components, which can bias their values. Our tests suggest indeed that better results are obtained (i.e. higher SNR for injected fake planets) when no flux scaling factor is considered before PCA.

7.5.3.1. Full-frame PCA

To trigger SDI alone with the psfsub.pca function, you need to set:

  • adimsdi='double' to separate the PCA-SDI and PCA-ADI
  • ncomp to a tuple of 2 values, the first of which being the number of principal components for PCA-SDI, and the second to be set to None - to prevent PCA-ADI (i.e. simply derotate and stack the SDI residual images);

Let’s use the mask_center_px argument to set a numerical mask of 0.1’’ radius in order not to include the coronagraphic mask.

[40]:

ncomp_sdi = 1
mask_px = int(0.1/plsc)

pca_sdi = pca(cube_fc, derot_angles, scale_list=opt_scal_med, ncomp=(ncomp_sdi,None),
              adimsdi='double', crop_ifs=False, imlib='opencv',
              interp_zeros=True, mask_val=0, mask_center_px=mask_px)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:48:28
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 0.938 GB
39 spectral channels in IFS cube
First PCA stage exploiting spectral variability
Running time:  0:00:16.624008
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
De-rotating and combining frames (skipping PCA)
Running time:  0:00:16.676963
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s show the final image with and without the application of a large (0.15’’ radius) numerical mask, and with circles at the ground truth location where we injected the fake companions.

[41]:

%matplotlib inline
plot_frames((pca_sdi, mask_circle(pca_sdi,20)),
            circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_105_0.png

In the first image, we notice the same effect as with median-SDI: significant stellar residuals at the edge of the coronagraphic mask.

On the other hand, fake planet c is easily re-detected. Notice the negative side lobe, that is reminiscent of ADI side lobes.

Question 7.1: Why is the negative side lobe only radially inward of the planet in our final SDI image?

This time the images obtained without and with flux-scaling applied to the cube before PCA appears to yield similar results. We note that flux scaling may actually have been slightly counter-productive, based on the raw fluxes at the location of the companion (peak value closer to ~0.4 without flux scaling, while closer to ~0.3 with flux scaling). Let’s confirm that by measuring the SNR of c:

[42]:

snr_c = snr(pca_sdi, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c)

13.017909840326496

7.5.3.2. Annular PCA

PCA-SDI can also be performed in concentric annuli using psfsub.pca_annular. Again let’s demonstrate whether flux scaling is useful or not:

[43]:

ncomp_sdi = 1
mask_px = int(0.1/plsc)
fwhm_m = np.mean(fwhm)

pca_sdi_ann = pca_annular(cube_fc, derot_angles, scale_list=opt_scal_med,
                          ncomp=(ncomp_sdi,None), radius_int=mask_px, fwhm=fwhm_m, asize=fwhm_m,
                          delta_sep=(0.4, 1), delta_rot=(0.5, 1),
                          mask_val=0, interp_zeros=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:48:45
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
First PCA subtraction exploiting the spectral variability
39 spectral channels per IFS frame
N annuli = 6, mean FWHM = 6.000
Running time:  0:00:33.878158
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
Skipping the second PCA subtraction

We mask the central part of the image, dominated by residual stellar light, to only focus on planet c:

[44]:

plot_frames(mask_circle(pca_sdi_ann, int(0.18/plsc)), circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_114_0.png 
[45]:

snr_c = snr(pca_sdi_ann, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c)

16.07169997049233

Despite both reductions yielding relatively similar final images, we note that flux-scaling is slightly counter-productive (SNR slightly lower), similarly to the full-frame PCA reduction.

7.5.4. median-ASDI

Let’s now combine both SDI and ADI, which should yield a better subtraction of the stellar residuals. Again, let’s start with the most basic way to combine both strategies: median-ASDI.

This simply involves the same calls to median_sub as in Sec. 7.5.2., but removing sdi_only=True:

[46]:

mask_px = int(0.1/plsc)

med_ASDI = median_sub(cube_fc, derot_angles, scale_list=opt_scal_med, flux_sc_list=opt_flux_med,
                      radius_int=mask_px, interp_zeros=True, mask_val=0)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:49:21
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
39 spectral channels per IFS frame
First median subtraction exploiting spectral variability
Running time:  0:00:17.307584
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
Median subtraction in the ADI fashion
Done derotating and combining
Running time:  0:00:19.472696
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s check the ASDI image:

[47]:

%matplotlib inline
plot_frames(med_ASDI, circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_121_0.png

This time we find signal at the location of fake planet b - although it does appear significantly filtered.

Fake planet c is easily redetected. Let’s measure the signal-to-noise ratio for both fake planets:

[48]:

snr_b = snr(med_SDI, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_b)

snr_c = snr(med_SDI, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c)

0.9083706649570626
19.808257154238373

7.5.5. PCA-ASDI

7.5.4.1. Full-frame PCA-ASDI (single step)

The psfsub.pca function in VIP offers the possibility to do this in either a single or 2 consecutive steps.

In the first case, a master PCA library is created with all images acquires at different parallactic angles and at different wavelengths (i.e. containing both angular and radial diversity). This is done by setting:

  • adimsdi='single' (default);
  • ncomp to a scalar value (integer to use that number of principal components, or float betweeb 0 and 1 to indicate the desired cumulative explained variance ratio).
[49]:

ncomp = 6
pca_asdi = pca(cube_fc, derot_angles, scale_list=opt_scal_med, ncomp=ncomp,
                adimsdi='single', crop_ifs=False, mask_center_px=mask_px,
                interp_zeros=True, mask_val=0)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:49:41
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 0.852 GB
Rescaling the spectral channels to align the speckles
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:00:08
Running time:  0:00:08.329935
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2145 total frames
Performing single-pass PCA
Done vectorizing the frames. Matrix shape: (2145, 22201)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:18.353620
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Descaling the spectral channels
0% [##############################] 100% | ETA: 00:00:00
Total time elapsed: 00:00:06
De-rotating and combining residuals
Running time:  0:00:24.546848
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Let’s visualize the resulting image, including circles at the ground truth locations of the injected planets:

[50]:

plot_frames(pca_asdi, label='PCA-ASDI (npc={:.0f})'.format(ncomp),
            dpi=100, colorbar=True, circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_129_0.png

Let’s estimate the signal-to-noise ratio of both detections:

[51]:

snr_b = snr(pca_asdi, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_b)

snr_c = snr(pca_asdi, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c)

5.622033970466388
27.689109243199827

Again we notice that pre-scaling the cube before PCA does not help with the detections.

7.5.4.2. Full-frame PCA-ASDI (double step)

The second option consists in performing PCA-SDI first, and then PCA-ADI on the cube consisting of SDI frames obtained for each individual spectral cube. This mode can be activated by setting adimsdi='double' (instead of the default ‘single’ value):

[52]:

ncomp = (1,1)
pca_asdi_2s = pca(cube_fc, derot_angles, scale_list=opt_scal_med, ncomp=ncomp,
                adimsdi='double', crop_ifs=False, mask_center_px=mask_px,
                interp_zeros=True, mask_val=0)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:50:08
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 1.070 GB
39 spectral channels in IFS cube
First PCA stage exploiting spectral variability
Running time:  0:00:16.261237
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
Second PCA stage exploiting rotational variability
De-rotating and combining residuals
Running time:  0:00:18.420811
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[53]:

plot_frames(pca_asdi_2s,
            label='PCA-ASDI 2s., npc=({:.0f},{})'.format(ncomp[0],ncomp[1]),
            dpi=100, colorbar=True, circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_136_0.png

Let’s estimate the signal-to-noise ratio of both detections:

[54]:

snr_b = snr(pca_asdi_2s, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_b)

snr_c = snr(pca_asdi_2s, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c)

5.078851917199891
50.37445233216146

7.5.4.3. Annular PCA-ASDI (double step)

[55]:

ncomp = (1,1)
pca_ann = pca_annular(cube_fc, derot_angles, scale_list=opt_scal_med, ncomp=ncomp,
                      radius_int=mask_px, fwhm=np.mean(fwhm), asize=np.mean(fwhm),
                      mask_val=0, interp_zeros=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:50:27
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
First PCA subtraction exploiting the spectral variability
39 spectral channels per IFS frame
N annuli = 6, mean FWHM = 6.000
Running time:  0:00:31.752404
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
55 ADI frames
Second PCA subtraction exploiting the angular variability
N annuli = 6, FWHM = 6.000
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  2.15    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:31.854678
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  4.36    Ann center:  22    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:31.927668
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  5.63    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:32.019417
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 4    PA thresh:  6.44    Ann center:  34    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:32.076253
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 5    PA thresh:  7.01    Ann center:  40    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:32.195415
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 6    PA thresh:  7.60    Ann center:  45    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:32.325645
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:34.451524
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[56]:

plot_frames(pca_ann, label='PCA-ASDI ann. (npc={:.0f},{:.0f})'.format(ncomp[0], ncomp[1]),
            dpi=100, colorbar=True, circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_141_0.png

Let’s estimate the signal-to-noise ratio and signifance of both detections:

[57]:

snr_b = snr(pca_ann, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_b)

snr_c = snr(pca_ann, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c)

7.538455403272426
38.36549205878312

[58]:

sig_b = significance(snr_b, r_b, np.mean(fwhm))
print(sig_b)

sig_c = significance(snr_c, r_c, np.mean(fwhm))
print(sig_c)

4.775173036007135
Warning high S/N! cdf>0.9999999999999999 is rounded to 1
Returning 8.2 sigma, but quote significance > 8.2 sigma.
8.2

7.5.6. PCA-ADI for all spectral channels

If you provide a 4d input cube without scale_list to the pca or pca_annular functions, their default behaviour will be to perform PCA-ADI on each spectral channel and then combine the final ADI images from each spectral channel into a single one.

[59]:

mask_px = int(0.08/plsc)

Let’s try different flavours of PCA. First let’s run it in full frames:

[60]:

ncomp_ff = 6
pca_adi_mean = pca(cube_fc, derot_angles, ncomp=ncomp_ff, adimsdi='single', mask_center_px=mask_px,
                   interp_zeros=True, mask_val=0)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:51:01
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
System total memory = 17.180 GB
System available memory = 0.957 GB
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.023710
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:02.157702
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:02.186426
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:04.310738
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:04.335221
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:06.468478
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:06.505410
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:08.653031
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:08.673518
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:10.800946
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:10.860721
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:12.993896
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:13.026665
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:15.209738
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:15.235866
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:17.370091
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:17.396284
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:19.521641
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:19.552033
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:21.676474
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:21.707453
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:23.829039
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:23.861539
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:25.989594
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:26.010170
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:28.158281
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:28.193291
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:30.344954
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:30.373778
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:32.524232
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:32.564310
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:34.700978
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:34.722973
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:36.863127
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:36.896668
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:39.034029
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:39.063261
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:41.205472
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:41.235773
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:43.373442
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:43.392773
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:45.528239
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:45.554470
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:47.689596
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:47.714318
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:49.851389
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:49.874103
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:52.022748
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:52.051989
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:54.190102
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:54.215494
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:56.364962
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:56.392223
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:00:58.526108
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:58.545358
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:00.702668
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:00.724572
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:02.880153
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:02.899594
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:05.040907
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:05.058395
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:07.232503
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:07.256528
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:09.429351
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:09.514646
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:11.668822
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:11.711320
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:13.866213
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:14.024321
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:16.209349
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:16.232011
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:18.395506
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:18.456014
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:20.618518
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:20.638436
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:22.788819
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done vectorizing the frames. Matrix shape: (55, 11449)
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:01:22.821789
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done de-rotating and combining
Running time:  0:01:24.986152
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Now in concentric annuli, with a PA threshold (set by delta_rot):

[61]:

ncomp_ann = 4
pca_adi_ann_mean = pca_annular(cube_fc, derot_angles, ncomp=ncomp_ann, radius_int=mask_px,
                               fwhm=np.mean(fwhm), asize=2*np.mean(fwhm), delta_rot=(0.2,1),
                               mask_val=0, interp_zeros=True)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:52:26
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.145423
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.401000
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:00.693143
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:02.874146
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:02.963390
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:03.341022
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:03.787270
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:05.931113
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:06.018148
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:06.162225
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:06.361933
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:08.497122
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:08.587040
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:08.756354
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:08.971969
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:11.121326
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:11.191089
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:11.369597
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:11.558599
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:13.695115
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:13.783205
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:13.963715
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:14.140011
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:16.292764
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:16.393253
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:16.563325
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:16.732668
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:18.875979
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:18.965474
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:19.117925
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:19.300678
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:21.442088
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:21.516755
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:21.652392
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:21.820470
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:23.947682
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:24.035248
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:24.199867
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:24.396152
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:26.520712
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:26.591773
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:26.756245
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:26.902113
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:29.021554
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:29.095604
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:29.253398
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:29.403710
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:31.550530
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:31.634915
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:31.853616
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:32.014749
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:34.138289
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:34.233911
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:34.407332
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:34.547381
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:36.670621
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:36.744259
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:36.873607
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:37.030026
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:39.154255
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:39.261335
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:39.428116
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:39.593911
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:41.755568
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:41.844846
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:42.026909
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:42.225104
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:44.372261
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:44.475560
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:44.651510
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:44.854321
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:46.992052
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:47.089030
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:47.260491
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:47.458526
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:49.597716
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:49.684912
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:49.864334
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:50.053340
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:52.207509
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:52.279357
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:52.429722
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:52.569696
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:54.709878
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:54.793174
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:54.954602
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:55.108311
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:57.230793
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:57.309441
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:57.469003
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:57.618614
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:00:59.739758
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:59.827258
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:00:59.986928
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:00.129737
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:02.252210
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:02.336948
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:02.485081
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:02.627315
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:04.756843
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:04.837282
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:05.041944
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:05.217024
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:07.352268
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:07.432114
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:07.617594
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:07.845759
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:09.980058
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:10.067426
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:10.301605
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:10.536422
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:12.664975
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:12.744911
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:12.900497
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:13.124614
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:15.273669
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:15.374349
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:15.539124
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:15.712355
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:17.866127
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:17.973720
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:18.186177
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:18.344339
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:20.485211
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:20.566368
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:20.780040
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:20.961436
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:23.100068
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:23.232044
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:23.491619
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:23.705945
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:25.859082
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:25.968951
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:26.155522
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:26.365229
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:28.510020
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:28.586404
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:28.783601
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:28.931489
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:31.063925
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:31.139801
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:31.291809
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:31.469575
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:33.597009
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:33.662794
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:33.825194
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:34.036409
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:36.167991
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:36.254524
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:36.410166
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:36.568145
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:38.693294
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
N annuli = 3, FWHM = 6.031
PCA per annulus (or annular sectors):
Ann 1    PA thresh:  4.31    Ann center:  16    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:38.775202
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 2    PA thresh:  7.37    Ann center:  28    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:38.929791
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Ann 3    PA thresh:  8.81    Ann center:  39    N segments: 1
Done PCA with lapack for current annulus
Running time:  0:01:39.090887
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done derotating and combining.
Running time:  0:01:41.218660
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Finally, let’s also run it as a single wide annulus including both injected planets:

[62]:

ncomp_a = 3
pca_adi_a_mean = pca_annulus(cube_fc, derot_angles, ncomp=ncomp_a, r_guess=(r_b+r_c)/2,
                           annulus_width=r_c, mask_val=0, interp_zeros=True)

Let’s compare the results:

[63]:

plot_frames((pca_adi_mean, pca_adi_ann_mean, pca_adi_a_mean),
            label=('PCA-ADI (npc={:.0f})'.format(ncomp_ff),
                   'PCA-ADI annular (npc={:.0f})'.format(ncomp_ann),
                   'PCA-ADI annulus (npc={:.0f})'.format(ncomp_a)),
            dpi=100, colorbar=True, circle=(xy_b,xy_c))
../_images/tutorials_07_ifs_psfsub_fm_planets_155_0.png

Let’s estimate the signal-to-noise ratio and signifance of both detections:

[64]:

snr_b_ff = snr(pca_adi_mean, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
snr_b_ann = snr(pca_adi_ann_mean, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
snr_b_a = snr(pca_adi_a_mean, xy_b, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_b_ff, snr_b_ann, snr_b_a)
snr_b = max(snr_b_ff, snr_b_ann, snr_b_a)

snr_c_ff = snr(pca_adi_mean, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
snr_c_ann = snr(pca_adi_ann_mean, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
snr_c_a = snr(pca_adi_a_mean, xy_c, fwhm=np.mean(fwhm), exclude_negative_lobes=True)
print(snr_c_ff, snr_c_ann, snr_c_a)
snr_c = max(snr_c_ff, snr_c_ann, snr_c_a)

9.93427387990257 6.2099793411402615 4.704562950646762
31.480749701607085 11.3982404029109 24.704716561264384

[65]:

sig_b = significance(snr_b, r_b, np.mean(fwhm))
print(sig_b)

sig_c = significance(snr_c, r_c, np.mean(fwhm))
print(sig_c)

5.433136382174582
Warning high S/N! cdf>0.9999999999999999 is rounded to 1
Returning 8.2 sigma, but quote significance > 8.2 sigma.
8.2

In this case, we see that ADI alone recovers better the innermost planet. This is because of the bright stellar residuals at the edge of the coronagraphic mask, preventing SDI to do a good job at good separation.

Answer 7.1: The scaling factors are such that all spectral channels are upscaled, for the stellar halo to match the last spectral channel, i.e. all factors are >= 1 (exactly 1 for the last spectral channel). The PSF model and subtraction with PCA is then made with these upscaled images, where the planet is now present at different radial separations. However since we injected a very red spectrum for planet c (blackbody of 1100 K), it is the brightest in the last spectral channels. This

means the strongest positive (resp. negative) imprint will be in the long (resp. short) wavelength channels, where the planet was not moved (resp. moved radially outward), i.e. the positive (resp. negative) imprint is unchanged (resp. radially inward) with respect to its position in the upscaled cube. After final downscaling, the negative imprint stays thus radially inward from the location of the fake planet. You can test this explanation by injecting a blue spectrum for planet c (e.g. 1/blackbody).

7.6. Astrometry retrieval

The firstguess function can accept a 4D input cube to infer a first guess on the position (a single one for all channels) and fluxes (one per channel). Since it would be very slow to run it on all 39 spectral channels, let’s just focus on the 8 channels (~20% of all channels) with the maximum S/N ratio. Let’s find these channels for planet c.

As the rest of this tutorial is CPU-intensive, we will only focus on planet c - although the same procedure can be applied to planet b.

7.6.1. Spectral channels with maximum SNR

[66]:

all_frs = []
all_snrs = []
all_pcs = []
for i in range(nch):
    res_ann_opt = pca_grid(cube_fc[i], derot_angles, fwhm=fwhm[i], range_pcs=(1,11,1),
                           source_xy=xy_c, mode='annular', annulus_width=int(4*fwhm[i]),
                           imlib='opencv', full_output=True, plot=True, exclude_negative_lobes=True)
    _, final_ann_opt, df, opt_npc_ann = res_ann_opt
    all_frs.append(final_ann_opt)
    pcs = df['PCs']
    snrs = df['S/Ns']
    opt_snr = snrs[df.index[pcs==opt_npc_ann]]
    print("Opt SNR: {}, for npc={}".format(float(opt_snr), opt_npc_ann))
    all_snrs.append(float(opt_snr))
    all_pcs.append(opt_npc_ann)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:33
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.028333
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=1.380
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 0.213
Central pixel S/N = 1.124
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 1.375
Max S/N (shifting the aperture center) = 2.200
stddev S/N (shifting the aperture center) = 0.565

Opt SNR: 1.3802059637391422, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:34
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.026543
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=3.570
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.134
Central pixel S/N = 4.537
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 3.582
Max S/N (shifting the aperture center) = 4.669
stddev S/N (shifting the aperture center) = 0.761

Opt SNR: 3.5696273024187377, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:36
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.026909
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=3.384
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.378
Central pixel S/N = 3.798
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 3.380
Max S/N (shifting the aperture center) = 4.865
stddev S/N (shifting the aperture center) = 0.671

Opt SNR: 3.3844253986246327, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:37
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.052070
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=3.538
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.718
Central pixel S/N = 4.374
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 3.554
Max S/N (shifting the aperture center) = 4.297
stddev S/N (shifting the aperture center) = 0.503

Opt SNR: 3.5377664065477767, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:38
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.026627
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=3.486
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.989
Central pixel S/N = 4.425
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 3.488
Max S/N (shifting the aperture center) = 4.297
stddev S/N (shifting the aperture center) = 0.507

Opt SNR: 3.48607005725874, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:39
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.026728
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=3.913
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.236
Central pixel S/N = 4.852
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 3.906
Max S/N (shifting the aperture center) = 4.805
stddev S/N (shifting the aperture center) = 0.594

Opt SNR: 3.9129698309141645, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.038780
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=3.632
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.046
Central pixel S/N = 4.416
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 3.601
Max S/N (shifting the aperture center) = 4.576
stddev S/N (shifting the aperture center) = 0.568

Opt SNR: 3.6323420173702496, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:41
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.041293
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=4.143
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.449
Central pixel S/N = 5.386
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.148
Max S/N (shifting the aperture center) = 5.303
stddev S/N (shifting the aperture center) = 0.670

Opt SNR: 4.143193146064415, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:42
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.030953
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=4.457
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.664
Central pixel S/N = 5.837
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.441
Max S/N (shifting the aperture center) = 5.895
stddev S/N (shifting the aperture center) = 0.880

Opt SNR: 4.457369370112796, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:43
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.013950
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=4.096
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.551
Central pixel S/N = 5.117
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.025
Max S/N (shifting the aperture center) = 5.454
stddev S/N (shifting the aperture center) = 0.851

Opt SNR: 4.095777478315322, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:44
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.027842
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=4.282
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 3.074
Central pixel S/N = 5.362
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.252
Max S/N (shifting the aperture center) = 5.566
stddev S/N (shifting the aperture center) = 0.937

Opt SNR: 4.2817197759528005, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:45
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.027366
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=4.445
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 0.997
Central pixel S/N = 5.858
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.115
Max S/N (shifting the aperture center) = 6.379
stddev S/N (shifting the aperture center) = 1.112

Opt SNR: 4.444660506253272, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:46
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.021464
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=5.323
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 0.983
Central pixel S/N = 6.558
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.820
Max S/N (shifting the aperture center) = 7.051
stddev S/N (shifting the aperture center) = 1.409

Opt SNR: 5.3233022967167365, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:47
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.127306
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=6.188
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.252
Central pixel S/N = 7.141
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 5.329
Max S/N (shifting the aperture center) = 7.451
stddev S/N (shifting the aperture center) = 1.474

Opt SNR: 6.1878561657941775, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:49
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.019968
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 6, for S/N=6.063
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.062
Central pixel S/N = 7.120
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 5.028
Max S/N (shifting the aperture center) = 7.533
stddev S/N (shifting the aperture center) = 1.633

Opt SNR: 6.063139538957518, for npc=6
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:50
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.017301
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 7, for S/N=7.297
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.233
Central pixel S/N = 8.846
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 6.375
Max S/N (shifting the aperture center) = 9.764
stddev S/N (shifting the aperture center) = 1.999

Opt SNR: 7.296945234614729, for npc=7
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:51
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.035918
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=5.542
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.950
Central pixel S/N = 7.021
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 5.521
Max S/N (shifting the aperture center) = 8.210
stddev S/N (shifting the aperture center) = 1.461

Opt SNR: 5.542413637737122, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:52
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.029605
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=4.976
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.305
Central pixel S/N = 5.791
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 5.041
Max S/N (shifting the aperture center) = 7.733
stddev S/N (shifting the aperture center) = 1.357

Opt SNR: 4.976404329522146, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:53
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.026733
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=4.888
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.248
Central pixel S/N = 6.085
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 4.908
Max S/N (shifting the aperture center) = 7.287
stddev S/N (shifting the aperture center) = 1.327

Opt SNR: 4.887619434458941, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:54
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.025238
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 2, for S/N=6.132
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.768
Central pixel S/N = 8.101
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 5.968
Max S/N (shifting the aperture center) = 8.959
stddev S/N (shifting the aperture center) = 1.641

Opt SNR: 6.13170533873802, for npc=2
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:55
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.018504
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 3, for S/N=7.050
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 1.945
Central pixel S/N = 8.241
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 6.651
Max S/N (shifting the aperture center) = 9.537
stddev S/N (shifting the aperture center) = 1.668

Opt SNR: 7.049993530235283, for npc=3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:56
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.015249
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

/Users/valentin/GitHub/VIP/vip_hci/psfsub/utils_pca.py:373: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.figure(figsize=vip_figsize, dpi=vip_figdpi)

Number of steps 11
Optimal number of PCs = 3, for S/N=8.520
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.583
Central pixel S/N = 9.774
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 7.872
Max S/N (shifting the aperture center) = 12.441
stddev S/N (shifting the aperture center) = 2.382

Opt SNR: 8.519994237380882, for npc=3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:57
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.020620
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=10.523
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.150
Central pixel S/N = 12.074
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 8.799
Max S/N (shifting the aperture center) = 12.194
stddev S/N (shifting the aperture center) = 2.371

Opt SNR: 10.522524996014564, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:58
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.025548
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=10.808
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.420
Central pixel S/N = 12.540
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 8.680
Max S/N (shifting the aperture center) = 12.377
stddev S/N (shifting the aperture center) = 2.430

Opt SNR: 10.807930290602918, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:55:59
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.094019
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=13.470
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.973
Central pixel S/N = 13.874
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.462
Max S/N (shifting the aperture center) = 14.511
stddev S/N (shifting the aperture center) = 2.743

Opt SNR: 13.469951242731133, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:01
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.024862
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=11.545
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.677
Central pixel S/N = 13.461
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 8.937
Max S/N (shifting the aperture center) = 13.917
stddev S/N (shifting the aperture center) = 3.184

Opt SNR: 11.545261112669461, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:02
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.021827
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=12.400
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 2.907
Central pixel S/N = 14.168
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 9.529
Max S/N (shifting the aperture center) = 14.563
stddev S/N (shifting the aperture center) = 3.334

Opt SNR: 12.399713012235242, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:03
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.017387
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=14.401
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 3.256
Central pixel S/N = 14.395
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.335
Max S/N (shifting the aperture center) = 16.271
stddev S/N (shifting the aperture center) = 3.354

Opt SNR: 14.401109970140446, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:04
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.019659
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=15.905
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 3.374
Central pixel S/N = 13.678
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.020
Max S/N (shifting the aperture center) = 14.898
stddev S/N (shifting the aperture center) = 3.339

Opt SNR: 15.905465591981226, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:05
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.024675
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=16.850
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 3.725
Central pixel S/N = 13.696
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.345
Max S/N (shifting the aperture center) = 15.549
stddev S/N (shifting the aperture center) = 3.404

Opt SNR: 16.85047458726297, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:06
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.023376
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=18.350
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 3.657
Central pixel S/N = 13.588
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.839
Max S/N (shifting the aperture center) = 16.412
stddev S/N (shifting the aperture center) = 3.782

Opt SNR: 18.34990761597011, for npc=5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:07
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.015516
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=17.783
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 4.466
Central pixel S/N = 17.642
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 12.035
Max S/N (shifting the aperture center) = 18.667
stddev S/N (shifting the aperture center) = 3.600

Opt SNR: 17.783037271804698, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:08
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.020120
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 3, for S/N=18.553
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 7.253
Central pixel S/N = 18.294
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 12.658
Max S/N (shifting the aperture center) = 18.712
stddev S/N (shifting the aperture center) = 3.802

Opt SNR: 18.55280700722969, for npc=3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:10
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.022752
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 3, for S/N=19.464
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 7.656
Central pixel S/N = 17.124
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 11.872
Max S/N (shifting the aperture center) = 18.229
stddev S/N (shifting the aperture center) = 3.599

Opt SNR: 19.46401148374845, for npc=3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:11
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.019907
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 4, for S/N=21.409
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 5.131
Central pixel S/N = 18.370
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 13.539
Max S/N (shifting the aperture center) = 22.127
stddev S/N (shifting the aperture center) = 4.682

Opt SNR: 21.408869859050355, for npc=4
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:12
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.022393
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 3, for S/N=20.230
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 7.500
Central pixel S/N = 16.489
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 11.942
Max S/N (shifting the aperture center) = 18.383
stddev S/N (shifting the aperture center) = 3.846

Opt SNR: 20.22969108248692, for npc=3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:13
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.023100
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=19.659
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 16.885
Central pixel S/N = 21.916
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 16.410
Max S/N (shifting the aperture center) = 22.823
stddev S/N (shifting the aperture center) = 4.146

Opt SNR: 19.6590426644281, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:14
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.018671
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 1, for S/N=20.985
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 15.664
Central pixel S/N = 24.092
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 17.044
Max S/N (shifting the aperture center) = 24.572
stddev S/N (shifting the aperture center) = 4.642

Opt SNR: 20.98510093582345, for npc=1
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:15
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.020661
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 2, for S/N=23.447
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 29.9
Flux in a centered 1xFWHM circular aperture = 7.526
Central pixel S/N = 14.151
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 14.491
Max S/N (shifting the aperture center) = 19.392
stddev S/N (shifting the aperture center) = 3.315

Opt SNR: 23.44685044744188, for npc=2
../_images/tutorials_07_ifs_psfsub_fm_planets_164_3.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_4.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_5.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_6.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_7.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_8.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_9.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_10.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_11.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_12.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_13.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_14.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_15.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_16.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_17.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_18.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_19.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_20.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_21.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_22.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_23.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_24.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_25.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_26.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_27.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_28.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_29.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_30.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_31.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_32.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_33.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_34.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_35.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_36.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_37.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_38.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_39.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_40.png
../_images/tutorials_07_ifs_psfsub_fm_planets_164_41.png

Let’s plot the results:

[67]:

plt.figure(figsize=(6,4))
plt.plot(lbda, all_snrs)
plt.xlabel(r"Wavelength ($\mu$m)")
plt.ylabel("SNR of c")

[67]:

Text(0, 0.5, 'SNR of c')
../_images/tutorials_07_ifs_psfsub_fm_planets_166_1.png

Let’s (reverse) sort to only keep the best 8 channels:

[68]:

order_snr = np.flip(np.argsort(all_snrs))
sorted_npc = np.array(all_pcs)[order_snr]
sorted_ch = np.arange(nch)[order_snr]
sorted_frs = np.array(all_frs)[order_snr]
sorted_snr = np.array(all_snrs)[order_snr]

Let’s finally show these 8 channels with optimal SNR:

[69]:

cutoff = 8
final_frs_c = tuple(sorted_frs[:cutoff])
final_chs_c = sorted_ch[:cutoff]
final_npcs_c = sorted_npc[:cutoff]

labels = ['ch.{}, npc={}, SNR={:.1f}'.format(final_chs_c[i], final_npcs_c[i], sorted_snr[i]) for i in range(cutoff)]
labels = tuple(labels)

plot_frames(final_frs_c, label=labels, rows=2)
../_images/tutorials_07_ifs_psfsub_fm_planets_170_0.png

7.6.2. NEGFC+PCA-ADI (simplex)

Let’s now find the optimal position of planet c by running the firstguess function on the 8 best spectral channels (corresponding to the highest SNR, as determined above):

[70]:

crop_cube_c = cube_fc[final_chs_c]
crop_cube_c.shape

[70]:

(8, 55, 107, 107)

[71]:

psfn_crop_c = psfn[final_chs_c]

[72]:

trans_crop_c = [trans[0]]
for i in final_chs_c:
    trans_crop_c.append(trans[i])
trans_crop_c = np.array(trans_crop_c)

[73]:

res = firstguess(crop_cube_c, derot_angles, psfn_crop_c, ncomp=list(final_npcs_c),
                 planets_xy_coord=[xy_c], fwhm=np.mean(fwhm), annulus_width=int(4*np.mean(fwhm)),
                 aperture_radius=2, imlib='opencv', transmission=trans_crop_c, mu_sigma=True,
                 force_rPA=False, simplex=True, plot=True, verbose=True)
r0c = res[0]
theta0c = res[1]
fluxes0c = res[2]

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 18:56:24
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [83.8,29.9], running ...
Processing spectral channel 0...
... optimal grid flux: 25.929 ($\chi^2_r$ = 2.2)
Processing spectral channel 1...
... optimal grid flux: 25.929 ($\chi^2_r$ = 2.3)
Processing spectral channel 2...
... optimal grid flux: 25.929 ($\chi^2_r$ = 1.9)
Processing spectral channel 3...
... optimal grid flux: 38.566 ($\chi^2_r$ = 2.2)
Processing spectral channel 4...
... optimal grid flux: 25.929 ($\chi^2_r$ = 1.9)
Processing spectral channel 5...
... optimal grid flux: 25.929 ($\chi^2_r$ = 2.2)
Processing spectral channel 6...
... optimal grid flux: 25.929 ($\chi^2_r$ = 2.2)
Processing spectral channel 7...
... optimal grid flux: 17.433 ($\chi^2_r$ = 2.4)
../_images/tutorials_07_ifs_psfsub_fm_planets_176_1.png 
Planet 0: preliminary position guess: (r, theta)=(38.5, 323.1)
Planet 0: preliminary flux guess: 25.9, 25.9, 25.9, 38.6, 25.9, 25.9, 25.9, 17.4
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 537, nfev: 944, chi2r: 0.15758758133832682
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f0, f1, f2, f3, f4, f5, f6, f7)=(25.360, 26.380, 29.330, 32.543, 31.531, 27.846, 26.891, 16.504) at
          (X,Y)=(83.75, 29.95)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:12:50.346442
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[74]:

res

[74]:

(array([38.43293289]),
 array([323.14460964]),
 array([[25.35980634, 26.38040389, 29.32954159, 32.54273905, 31.53075666,
         27.84553272, 26.89055049, 16.50429857]]))

Let’s compare the estimates to their respective ground truth values:

[75]:

print("r (gt): {:.2f} VS. r (estimated): {:.2f}".format(r_c, r0c[0]))
print("theta (gt): {:.2f} VS. r (estimated): {:.2f}".format(theta_c%360, theta0c[0]))
for i, ch in enumerate(final_chs_c):
    print("flux ch.{} (gt): {:.2f} VS. r (estimated): {:.2f}".format(ch, flux_c_scal[ch], fluxes0c[0,i]))

r (gt): 38.50 VS. r (estimated): 38.43
theta (gt): 323.13 VS. r (estimated): 323.14
flux ch.38 (gt): 24.66 VS. r (estimated): 25.36
flux ch.34 (gt): 29.42 VS. r (estimated): 26.38
flux ch.37 (gt): 28.25 VS. r (estimated): 29.33
flux ch.35 (gt): 30.51 VS. r (estimated): 32.54
flux ch.36 (gt): 30.09 VS. r (estimated): 31.53
flux ch.33 (gt): 27.34 VS. r (estimated): 27.85
flux ch.32 (gt): 25.13 VS. r (estimated): 26.89
flux ch.30 (gt): 22.10 VS. r (estimated): 16.50

7.6.3. NEGFC+PCA-ADI (MCMC)

A better, albeit slower way, to infer the astrometry of each planet with uncertainties is to use an MCMC approach. Our specific case can be solved with the mcmc_negfc_sampling, which can, since VIP v1.2.3, also deal with 4D cubes. We refer to tutorial 5 for more details on how to set up the parameters for the MCMC algorithm.

[76]:

import pickle

lab='c'
ann_width = 4*np.mean(fwhm)
aperture_radius=2
imlib_rot='opencv'
interpolation = 'lanczos4'
nwalkers, itermin, itermax = (100, 200, 500)
conv_test, ac_c, ac_count_thr, check_maxgap = ('ac', 50, 1, 50)

algo_params = {'algo': pca_annulus,
               'annulus_width': ann_width,
               'svd_mode': 'lapack',
               'imlib': imlib_rot,
               'interpolation': interpolation}
conv_params = {'conv_test': conv_test,
               'ac_c': ac_c,
               'ac_count_thr': ac_count_thr,
               'check_maxgap': check_maxgap}
mcmc_params = {'nwalkers': nwalkers,
               'niteration_min': itermin,
               'niteration_limit': itermax,
               'bounds': None,
               'sigma':'spe',
               'nproc': 2}
negfc_params = {'mu_sigma': True,
                'aperture_radius': aperture_radius}
obs_params = {}

Note that when the flux variation is known throughout the observation (e.g. through a background star or satellite spots) or if you wish to take into account the effect of airmass, it is possible to provide a weights parameter for the injection of the mean flux to be scaled at the moment of injection in each image to follow in (this is a new addition since Christiaens et al. 2021).

[77]:

weights=X_fac

[78]:

# input to MCMC
init = [r0c[0], theta0c[0]]
for i in range(cutoff):
    init.append(fluxes0c[0,i])
initial_state = np.array(init)
obs_params['psfn']= psfn[final_chs_c]
obs_params['fwhm']= np.mean(fwhm[final_chs_c]) # REPLACE AFTER DEALING WITH FWHM format
obs_params['transmission']= trans_crop_c
algo_params['ncomp'] = list(final_npcs_c)

Warning: The following box is very computer intensive. It can be skipped if needed.

[79]:

# MCMC
chain = mcmc_negfc_sampling(cube_fc[final_chs_c], derot_angles, **obs_params, **algo_params, **negfc_params,
                            initial_state=initial_state, **mcmc_params, **conv_params, weights=weights,
                            display=True, verbosity=2, save=False, output_dir='./')

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 19:09:15
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
        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.01 and 0.06 respectively.
Beginning emcee Ensemble sampler...
emcee Ensemble sampler successful

Start of the MCMC run ...
Step  |  Duration/step (sec)  |  Remaining Estimated Time (sec)
0               27.75770                        13851.09155
1               27.57119                        13730.45511
2               27.53320                        13684.00040
3               27.75947                        13768.69762
4               28.03294                        13876.30332
5               27.99338                        13828.72824
6               28.00774                        13807.81533
7               27.74790                        13651.96680
8               27.69214                        13596.84074
9               28.75231                        14088.63288
10              31.49064                        15398.92296
11              31.36916                        15308.15203
12              31.43031                        15306.56000
13              31.70301                        15407.66286
14              31.37198                        15215.40836
15              31.85860                        15419.56095
16              30.77905                        14866.28356
17              33.09605                        15952.29803
18              30.59164                        14714.57740
19              31.24397                        14997.10704
20              31.94149                        15299.97323
21              31.40366                        15010.94948
22              31.57124                        15059.48196
23              30.87921                        14698.50253
24              29.37687                        13954.01277
25              30.15498                        14293.46242
26              29.15031                        13788.09474
27              28.44059                        13423.95990
28              27.65868                        13027.24016
29              29.90354                        14054.66239
30              31.26557                        14663.55467
31              31.65902                        14816.42323
32              32.15261                        15015.27074
33              32.89504                        15329.08957
34              31.03462                        14431.09644
35              30.16904                        13998.43502
36              33.00778                        15282.60445
37              31.28254                        14452.53579
38              31.66442                        14597.29716
39              31.20322                        14353.48166
40              31.43530                        14428.80224
41              30.29527                        13875.23366
42              30.98146                        14158.52768
43              30.04717                        13701.50861
44              29.40969                        13381.41032
45              28.44601                        12914.49036
46              29.08780                        13176.77295
47              28.56146                        12909.77947
48              28.24334                        12737.74724
49              27.89304                        12551.86800
50              27.88959                        12522.42456

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_1.png 
51              32.43614                        14531.39206
52              29.95056                        13387.89898
53              30.46228                        13586.17554
54              29.64382                        13191.49812
55              29.83398                        13246.28668
56              29.57185                        13100.32999
57              30.70013                        13569.45879
58              30.03402                        13245.00326
59              30.33919                        13349.24404
60              29.92266                        13136.04818
61              30.41680                        13322.55971
62              30.50917                        13332.50554
63              30.26542                        13195.72443
64              29.50079                        12832.84365
65              30.46552                        13222.03785
66              29.41223                        12735.49602
67              29.85512                        12897.41270
68              29.79351                        12841.00453
69              29.93063                        12870.16918
70              29.23401                        12541.39072
71              29.52572                        12637.00730
72              29.95671                        12791.51389
73              29.60724                        12612.68339
74              29.40575                        12497.44503
75              30.38527                        12883.35533
76              30.21552                        12781.16496
77              29.83437                        12590.10330
78              29.90534                        12590.14982
79              29.43361                        12362.11788
80              30.01121                        12574.69531
81              30.30869                        12669.03200
82              30.26439                        12620.24938
83              29.57070                        12301.41120
84              29.64780                        12303.83617
85              29.59414                        12251.97272
86              30.19328                        12469.82629
87              30.00609                        12362.50908
88              30.64098                        12593.44114
89              29.52210                        12104.05936
90              28.86771                        11806.89257
91              152.96469                       62409.59270
92              55.55723                        22611.79220
93              863.52699                       350591.95713
94              754.85426                       305715.97692
95              29.38355                        11870.95501
96              29.80338                        12010.76295
97              29.63476                        11913.17392
98              29.41648                        11796.00928
99              29.47986                        11791.94480
100             29.86237                        11915.08483

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_3.png 
101             31.53162                        12549.58476
102             29.89731                        11869.23405
103             30.05866                        11903.22857
104             29.81645                        11777.49814
105             29.46966                        11611.04565
106             28.86931                        11345.63922
107             29.62273                        11612.11134
108             29.53402                        11547.80338
109             29.42287                        11474.91735
110             30.19960                        11747.64596
111             29.97464                        11630.15993
112             30.08042                        11641.12293
113             29.43142                        11360.52619
114             29.60791                        11399.04381
115             28.91886                        11104.84109
116             29.76230                        11398.95975
117             29.92868                        11432.75576
118             30.27510                        11534.81234
119             30.22483                        11485.43426
120             30.26643                        11470.97697
121             30.48112                        11521.86260
122             29.82104                        11242.53283
123             29.54129                        11107.52617
124             29.18201                        10943.25262
125             29.78163                        11138.32850
126             29.22987                        10902.74338
127             29.05712                        10809.24901
128             29.45885                        10929.23335
129             30.39060                        11244.52348
130             29.89205                        11030.16571
131             30.09308                        11074.25344
132             28.94219                        10621.78446
133             29.55284                        10816.34054
134             29.15118                        10640.18143
135             29.38967                        10697.83988
136             29.37990                        10664.90406
137             29.29459                        10604.64267
138             30.11574                        10871.78033
139             29.57093                        10645.53516
140             29.68222                        10655.91698
141             29.30666                        10491.78321
142             29.70616                        10605.09948
143             28.61274                        10186.13686
144             30.37365                        10782.64575
145             30.49625                        10795.67321
146             26.82673                        9469.83745
147             27.33696                        9622.61133
148             27.44144                        9631.94474
149             27.31726                        9561.04135
150             27.84419                        9717.62406

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_5.png 
151             29.03476                        10104.09718
152             27.75868                        9632.26022
153             26.31230                        9104.05545
154             26.82013                        9252.94554
155             27.27548                        9382.76340
156             27.31965                        9370.63995
157             27.27287                        9327.32257
158             26.75747                        9124.29727
159             27.30167                        9282.56780
160             26.31143                        8919.57477
161             26.46715                        8945.89771
162             26.76765                        9020.69906
163             27.80832                        9343.59686
164             27.20448                        9113.50247
165             27.24556                        9100.01704
166             27.26167                        9078.13644
167             27.20799                        9033.05135
168             27.25524                        9021.48444
169             26.81029                        8847.39702
170             26.72099                        8791.20505
171             27.24539                        8936.48858
172             26.37014                        8623.03513
173             27.17599                        8859.37144
174             27.73016                        9012.30233
175             26.73468                        8662.03794
176             27.29291                        8815.60961
177             27.26643                        8779.79143
178             27.21898                        8737.29322
179             26.92987                        8617.56000
180             27.22661                        8685.28731
181             27.35069                        8697.51942
182             27.16467                        8611.19944
183             26.09644                        8246.47662
184             25.72337                        8102.86249
185             26.73671                        8395.32631
186             26.94828                        8434.81039
187             26.62539                        8307.12043
188             26.27138                        8170.39856
189             27.07296                        8392.61884
190             27.16052                        8392.60192
191             27.12188                        8353.53904
192             26.55909                        8153.63971
193             26.71706                        8175.42158
194             27.38515                        8352.46983
195             27.04958                        8223.07323
196             27.01961                        8186.94092
197             27.35877                        8262.34884
198             27.82145                        8374.25495
199             26.73504                        8020.51080
200             27.21448                        8137.12982

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_7.png 
Auto-corr tau/N = [0.07931158 0.09041509 0.0817418  0.10018691 0.08593508 0.08128431
 0.08215103 0.0859475  0.08465429 0.10082706]
tau/N <= 0.02 = [False False False False False False False False False False]

201             29.26284                        8720.32513
202             27.52704                        8175.52999
203             27.09341                        8019.64847
204             27.49052                        8109.70310
205             27.57460                        8106.93328
206             27.46305                        8046.67394
207             27.81699                        8122.56108
208             27.70843                        8063.15284
209             27.67004                        8024.31218
210             27.08957                        7828.88660
211             28.09803                        8092.23322
212             27.62788                        7929.20127
213             27.42908                        7844.71602
214             26.79417                        7636.33845
215             27.15911                        7713.18809
216             26.77881                        7578.40266
217             27.46606                        7745.42977
218             27.65401                        7770.77737
219             27.20735                        7618.05688
220             27.21664                        7593.44172
221             27.55823                        7661.18877
222             28.09057                        7781.08817
223             27.01840                        7457.07812
224             27.03565                        7434.80403
225             27.84452                        7629.39930
226             27.90345                        7617.64185
227             28.49287                        7750.05982
228             27.31811                        7403.20862
229             26.72527                        7215.82263
230             27.25414                        7331.36232
231             26.85153                        7196.20977
232             26.67830                        7123.10557
233             26.87938                        7149.91402
234             27.42665                        7268.06119
235             27.36533                        7224.44686
236             26.70842                        7024.31446
237             26.98830                        7070.93329
238             26.67461                        6962.07399
239             26.94507                        7005.71820
240             27.18125                        7039.94453
241             27.49097                        7092.67103
242             27.55024                        7080.41117
243             27.30764                        6990.75610
244             26.98029                        6879.97446
245             26.94951                        6845.17554
246             27.68526                        7004.37179
247             27.00626                        6805.57702
248             26.92158                        6757.31558
249             26.82203                        6705.50650
250             26.83527                        6681.98323

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_9.png 
Auto-corr tau/N = [0.08020549 0.08832579 0.07607597 0.10064228 0.08229707 0.0765978
 0.07899672 0.08030464 0.08190058 0.10342003]
tau/N <= 0.02 = [False False False False False False False False False False]

251             29.71557                        7369.46012
252             26.89422                        6642.87333
253             26.28921                        6467.14640
254             26.53583                        6501.27786
255             26.55232                        6478.76510
256             26.08790                        6339.35921
257             26.30041                        6364.70019
258             26.18536                        6310.67224
259             26.52084                        6365.00208
260             26.34281                        6295.93255
261             26.23519                        6243.97474
262             27.13520                        6431.04335
263             26.85824                        6338.54464
264             26.36574                        6195.94890
265             26.38790                        6174.76766
266             26.59788                        6197.30651
267             26.53625                        6156.41023
268             27.40761                        6331.15675
269             26.45969                        6085.72962
270             26.35750                        6035.86636
271             26.54435                        6052.11248
272             26.89415                        6104.97296
273             27.21139                        6149.77482
274             26.79811                        6029.57565
275             25.89633                        5800.77770
276             27.53375                        6140.02692
277             26.49965                        5882.92230
278             26.34610                        5822.48898
279             26.61589                        5855.49602
280             25.54397                        5594.13009
281             26.76575                        5834.93285
282             25.96851                        5635.16580
283             25.69314                        5549.71910
284             26.12242                        5616.32137
285             26.45189                        5660.70446
286             26.80026                        5708.45623
287             26.02966                        5518.28792
288             24.83487                        5240.15820
289             27.21115                        5714.34192
290             26.80686                        5602.63332
291             26.42061                        5495.48709
292             26.63082                        5512.57974
293             25.66736                        5287.47595
294             26.30762                        5393.06189
295             26.15802                        5336.23547
296             26.78285                        5436.91916
297             26.31158                        5314.93876
298             26.38083                        5302.54663
299             26.52480                        5304.96040
300             27.59879                        5492.15941

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_11.png 
Auto-corr tau/N = [0.08209354 0.08617746 0.07863792 0.09758879 0.07731056 0.07473813
 0.07579321 0.08080995 0.07902462 0.10438969]
tau/N <= 0.02 = [False False False False False False False False False False]

301             27.41745                        5428.65431
302             26.52913                        5226.23940
303             26.39359                        5173.14384
304             25.67717                        5007.04835
305             24.75539                        4802.54585
306             26.50844                        5116.12853
307             26.71039                        5128.39430
308             25.10753                        4795.53899
309             27.21447                        5170.74911
310             25.46513                        4812.90919
311             26.68042                        5015.91952
312             25.94567                        4851.84104
313             26.02854                        4841.30807
314             26.10274                        4829.00690
315             26.38028                        4853.97207
316             26.88335                        4919.65305
317             25.51508                        4643.74383
318             26.41527                        4781.16297
319             26.88177                        4838.71860
320             25.77347                        4613.45149
321             25.90371                        4610.86002
322             26.57412                        4703.62012
323             26.37715                        4642.37770
324             26.61967                        4658.44225
325             27.25491                        4742.35434
326             27.40088                        4740.35293
327             26.20612                        4507.45298
328             27.47348                        4697.96457
329             26.58192                        4518.92657
330             27.42452                        4634.74405
331             27.58790                        4634.76720
332             27.21630                        4545.12243
333             27.71245                        4600.26753
334             27.50046                        4537.57541
335             26.20650                        4297.86649
336             26.57753                        4332.13674
337             25.99504                        4211.19697
338             25.53307                        4110.82347
339             26.48638                        4237.82032
340             27.29270                        4339.53962
341             25.79820                        4076.11607
342             26.26874                        4124.19297
343             26.99813                        4211.70750
344             27.03986                        4191.17768
345             26.52657                        4085.09147
346             26.31263                        4025.83208
347             26.17161                        3978.08426
348             26.76722                        4041.85022
349             26.45523                        3968.28465
350             26.30793                        3919.88083

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_13.png 
Auto-corr tau/N = [0.08474313 0.08508626 0.07534517 0.09432212 0.07509943 0.07298475
 0.07203059 0.08107009 0.07620588 0.10438039]
tau/N <= 0.02 = [False False False False False False False False False False]

351             28.11183                        4160.55158
352             26.54087                        3901.50804
353             26.59776                        3883.27281
354             26.46921                        3838.03617
355             26.73382                        3849.66950
356             27.47991                        3929.62756
357             26.85989                        3814.10481
358             26.66928                        3760.36820
359             27.43996                        3841.59468
360             26.07638                        3624.61724
361             27.12042                        3742.61782
362             26.90831                        3686.43902
363             27.33340                        3717.34308
364             27.46589                        3707.89556
365             27.43453                        3676.22675
366             26.95601                        3585.14906
367             26.70196                        3524.65846
368             27.15323                        3557.07313
369             26.05063                        3386.58203
370             25.07513                        3234.69164
371             26.04355                        3333.57440
372             25.92749                        3292.79136
373             26.48846                        3337.54646
374             26.04619                        3255.77350
375             27.31642                        3387.23620
376             27.13692                        3337.84079
377             26.70471                        3257.97462
378             26.50267                        3206.82319
379             26.50625                        3180.75060
380             27.57346                        3281.24150
381             26.72258                        3153.26409
382             27.03958                        3163.63039
383             26.65644                        3092.14727
384             26.92054                        3095.86164
385             26.58228                        3030.37969
386             25.95398                        2932.79985
387             27.67545                        3099.65062
388             25.67350                        2849.75839
389             27.14431                        2985.87421
390             26.66325                        2906.29447
391             27.32825                        2951.45111
392             26.10364                        2793.08916
393             25.12613                        2663.36936
394             26.23379                        2754.54795
395             26.65739                        2772.36898
396             26.65839                        2745.81468
397             27.40636                        2795.44821
398             26.65395                        2692.04844
399             27.14764                        2714.76360
400             26.79788                        2652.98972

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_15.png 
Auto-corr tau/N = [0.08275684 0.07983403 0.07388097 0.08987802 0.07285759 0.06947194
 0.07137076 0.0789889  0.07411971 0.10180646]
tau/N <= 0.02 = [False False False False False False False False False False]

401             28.35803                        2779.08665
402             28.16204                        2731.71798
403             25.18101                        2417.37658
404             26.17816                        2486.92482
405             27.19539                        2556.36628
406             25.59635                        2380.46074
407             26.10626                        2401.77601
408             26.04313                        2369.92447
409             26.95836                        2426.25195
410             26.81892                        2386.88379
411             28.16307                        2478.34981
412             27.26402                        2371.96948
413             25.83276                        2221.61710
414             27.19368                        2311.46288
415             26.66222                        2239.62673
416             26.54604                        2203.32124
417             27.26639                        2235.84365
418             26.17825                        2120.43817
419             27.17274                        2173.81896
420             26.95540                        2129.47628
421             27.34017                        2132.53342
422             27.34579                        2105.62598
423             26.64437                        2024.97204
424             26.49606                        1987.20480
425             25.49989                        1886.99179
426             26.18812                        1911.73305
427             27.21800                        1959.69578
428             26.46241                        1878.83090
429             26.46278                        1852.39432
430             28.45498                        1963.39383
431             27.34936                        1859.75682
432             27.87234                        1867.44705
433             27.08684                        1787.73170
434             27.13254                        1763.61536
435             26.04158                        1666.66080
436             26.61778                        1676.92001
437             26.61649                        1650.22263
438             26.17904                        1596.92156
439             26.59122                        1595.47302
440             27.45465                        1619.82423
441             25.69078                        1490.06524
442             26.42673                        1506.32384
443             26.82917                        1502.43369
444             28.04229                        1542.32617
445             25.63912                        1384.51253
446             26.68383                        1414.24304
447             27.04749                        1406.46969
448             27.10587                        1382.39927
449             26.23308                        1311.65405
450             26.77895                        1312.16860

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_187_17.png 
Auto-corr tau/N = [0.07921997 0.07633829 0.07374483 0.08831227 0.07115975 0.06782997
 0.07082012 0.07662235 0.07262194 0.09720397]
tau/N <= 0.02 = [False False False False False False False False False False]

451             28.47420                        1366.76141
452             26.91317                        1264.91913
453             26.75994                        1230.95729
454             26.91740                        1211.28309
455             26.10325                        1148.54282
456             27.42841                        1179.42154
457             27.29005                        1146.18193
458             27.10734                        1111.40078
459             26.53834                        1061.53344
460             26.36182                        1028.11094
461             28.30804                        1075.70533
462             26.38449                        976.22606
463             27.50631                        990.22716
464             27.37324                        958.06322
465             26.82272                        911.97258
466             26.86066                        886.40161
467             25.48422                        815.49498
468             27.98256                        867.45924
469             25.72727                        771.81798
470             26.08131                        756.35808
471             26.59919                        744.77726
472             26.90152                        726.34109
473             26.41181                        686.70696
474             26.47773                        661.94337
475             27.46916                        659.25977
476             26.91592                        619.06621
477             26.18504                        576.07084
478             26.37792                        553.93624
479             27.00388                        540.07754
480             27.05034                        513.95644
481             27.72610                        499.06975
482             25.54045                        434.18767
483             26.60975                        425.75592
484             26.75160                        401.27403
485             26.16497                        366.30965
486             25.97715                        337.70290
487             26.33219                        315.98623
488             26.29202                        289.21223
489             25.98549                        259.85488
490             25.87163                        232.84466
491             26.25694                        210.05549
492             26.53267                        185.72866
493             25.91881                        155.51285
494             26.70040                        133.50198
495             25.77507                        103.10028
496             26.80062                        80.40186
497             27.52763                        55.05527
498             27.20336                        27.20336
499             27.00046                        0.00000
We have reached the limit # of steps without convergence
Running time:  4:19:57.150176
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

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.

[80]:

write=False

if write:
    output = {'chain':chain}
    with open('../datasets/MCMC_results_pl{}_top{}ch'.format(lab, cutoff), 'wb') as fileSave:
        pickle.dump(output, fileSave)

[81]:

with open('../datasets/MCMC_results_pl{}_top{}ch'.format(lab, cutoff), 'rb') as fi:
    myPickler = pickle.Unpickler(fi)
    mcmc_result = myPickler.load()

chain = mcmc_result['chain']

Let’s visualize the chain after setting labels for each parameter:

[82]:

labels = [r"$r$", r"$\theta$"]
for i, ch in enumerate(final_chs_c):
    labels.append(r"$f$ (ch.{:.0f})".format(ch))
show_walk_plot(chain, labels=labels)
../_images/tutorials_07_ifs_psfsub_fm_planets_192_0.png

Let’s visualize the corner plot after burning the first 30% of the chain:

[83]:

burnin = 0.3
show_corner_plot(chain, burnin=burnin, labels=labels)
../_images/tutorials_07_ifs_psfsub_fm_planets_194_0.png

To calculate 1-sigma (68%) confidence intervals on each parameter, let’s first flatten the chain:

[84]:

npar = len(initial_state)
isamples_flat = chain[:, int(chain.shape[1]//(1/burnin)):, :].reshape((-1,npar))

Then use the confidence function.

Below, the first set of plots show the estimate, ground truth value and 68% confidence interval (in shaded area).

By setting ‘gaussian_fit’ to True, a second set of plots is also shown, where a Gaussian function is fit to each posterior distribution in order to infer mean and standard deviation for each parameter.

Note that the function can accept a ground truth value gt, which will be shown in the plots if provided.

[85]:

gt = [r_c, theta_c%360]
for i, ch in enumerate(final_chs_c):
    gt.append(flux_c_scal[ch])

mu, sigma = confidence(isamples_flat, cfd=68, gaussian_fit=True, verbose=False,
                       gt=gt, save=False, title=True, labels=labels)
../_images/tutorials_07_ifs_psfsub_fm_planets_198_0.png

Overall we see that the estimated parameters are consistent with the ground truth values (shown with dashed blue lines).

[86]:

r_est = mu[0]
theta_est = mu[1]

xc, yc = frame_center(cube_fc)
x_est = xc + r_est*np.cos(np.deg2rad(theta_est))
y_est = yc + r_est*np.sin(np.deg2rad(theta_est))

xy_est = (x_est, y_est)

7.7. Spectrum retrieval

7.7.1. NEGFC+PCA-ADI (simplex)

Let’s use the negative fake companion technique (combined with PCA-ADI) to estimate the flux of the companion in each spectral channel. Here, we will just use the firstguess function in VIP which finds a first estimate on a grid of flux values followed by a simplex. Here we set force_rPA=True to pin the position during the NEGFC procedure (i.e. only allowing the fluxes at different wavelengths as free parameters).

[87]:

est_fluxes_c = np.zeros(nch)
for i in range(nch):
    trans_i = np.array([trans[0],trans[i]])
    res = firstguess(cube_fc[i], derot_angles, psfn[i], ncomp=4, planets_xy_coord=[xy_est],
                     fwhm=fwhm[i], annulus_width=int(4*fwhm[i]), aperture_radius=2, imlib='opencv',
                     f_range=None, transmission=trans_i, mu_sigma=True, force_rPA=True, weights=X_fac,
                     simplex=True, plot=False, verbose=False)
    _, _, est_fluxes_c[i] = res

Note that we arbitrarily assigned the number of principal components to 5 and 2 for the PCA-ADI algorithm used by NEGFC to measure the fluxes of the inner and outer fake planet. A more rigorous approach would consist in finding the optimal number of principal components for each spectral channel first, as done in Tutorial 5 (Sec. 5.2).

Let’s now compare the retrieved fluxes with and without weights to the ground truth values used for injection:

[88]:

%matplotlib inline
fig, axes = plt.subplots(1,2, figsize = (14,5))
ax1, ax2 = axes
ax1.plot(lbda, est_fluxes_c, 'bo', label='Planet c: estimated spectrum')
ax1.plot(lbda, flux_c_scal, 'k', label='Planet c: injected spectrum')
ax1.set_xlabel("Wavelength")
ax1.set_ylabel("Flux (ADU/s)")
ax1.legend()

ax2.plot(lbda, est_fluxes_c-flux_c_scal, 'bo', label='Planet b: estimated - injected spectrum')
ax2.plot(lbda, [0]*len(lbda), 'k--')
ax2.plot(lbda, [np.mean(est_fluxes_c-flux_c_scal)]*len(lbda), 'b-', label='Mean residuals')
ax2.set_xlabel("Wavelength")
ax2.set_ylabel("Flux (ADU/s)")
ax2.legend()
plt.show()
../_images/tutorials_07_ifs_psfsub_fm_planets_207_0.png

This approach does not yield uncertainties. For the latter, it is recommended to run the mcmc_negfc_sampling function.

7.7.2. NEGFC+PCA-ADI (MCMC)

One can use the MCMC approach to have both estimates and uncertainties on the flux of directly imaged companions. Given that this approach is computationally expensive, we will only apply it to 3 spectral channels.

We refer to Tutorial 5 (Sec. 5.3.3.) for more details on setting the MCMC parameters.

[89]:

ann_width = 4*np.mean(fwhm)
aperture_radius=2
imlib_rot='opencv'
interpolation = 'lanczos4'
nwalkers, itermin, itermax = (100, 200, 500)
conv_test, ac_c, ac_count_thr, check_maxgap = ('ac', 50, 1, 50)

algo_params = {'algo': pca_annulus,
               'annulus_width': ann_width,
               'svd_mode': 'lapack',
               'imlib': imlib_rot,
               'interpolation': interpolation}
conv_params = {'conv_test': conv_test,
               'ac_c': ac_c,
               'ac_count_thr': ac_count_thr,
               'check_maxgap': check_maxgap}
mcmc_params = {'nwalkers': nwalkers,
               'niteration_min': itermin,
               'niteration_limit': itermax,
               'bounds': None,
               'sigma':'spe',
               'nproc': 2}
negfc_params = {'mu_sigma': True,
                'aperture_radius': aperture_radius}
obs_params = {}

Important note: the pixel intensities in calibrated SPHERE cubes coming out from the DRH pipeline (as this one) are normalized by their integration time. In absence of knowledge of this original integration time, it is not possible to scale the cube bask to original ADU values, and it is therefore not possible to include photon noise uncertainty among the different sources of uncertainty (in this specific case, the estimated photon noise uncertainty will be so large that the MCMC will not work properly). To alleviate this issue, we set sigma to ‘spe’, instead of ‘spe+pho’ (default) - although this may slightly underestimate the final uncertainties.

[90]:

import pickle

lab='c'
test_ch = [20,29,38]

Warning: the next box is computer-intensive. The results from that box have been saved in the ‘datasets’ folder, such that it can be skipped.

[91]:

for i, ch in enumerate(test_ch):
    # find optimal npc
    res_ann_opt = pca_grid(cube_fc[ch], derot_angles, fwhm=fwhm[ch], range_pcs=(1,11,1),
                           source_xy=xy_est, mode='annular',
                           annulus_width=ann_width, imlib=imlib_rot,
                           interpolation=interpolation,
                           full_output=True, plot=True, exclude_negative_lobes=True)
    _, final_ann_opt, _, opt_npc_ann = res_ann_opt
    plot_frames(final_ann_opt, label='PCA-ADI annulus (ch={}, opt. npc={:.0f})'.format(ch, opt_npc_ann),
                dpi=100, colorbar=True, circle=xy_est)
    # first guess
    trans_i = np.array([trans[0],trans[ch]])
    r_0, theta_0, f_0 = firstguess(cube_fc[ch], derot_angles, psfn[ch], ncomp=opt_npc_ann,
                                   planets_xy_coord=[xy_est], fwhm=fwhm[ch], weights=X_fac,
                                   f_range=None, annulus_width=ann_width, aperture_radius=aperture_radius,
                                   imlib=imlib_rot, interpolation=interpolation, simplex=True, force_rPA=True,
                                   transmission=trans_i, plot=True, verbose=True)
    initial_state = np.array([r_0[0], theta_0[0], f_0[0]])
    # MCMC
    obs_params['psfn']= psfn[ch]
    obs_params['fwhm']= fwhm[ch]
    obs_params['transmission']= trans_i
    algo_params['ncomp'] = opt_npc_ann
    chain = mcmc_negfc_sampling(cube_fc[ch], derot_angles, **obs_params, **algo_params, **negfc_params,
                                initial_state=initial_state, **mcmc_params, **conv_params, weights=X_fac,
                                force_rPA=True, display=True, verbosity=2, save=False, output_dir='./')
    output = {'chain':chain}
    with open('../datasets/MCMC_results_pl{}_ch{}'.format(lab, ch), 'wb') as fileSave:
        pickle.dump(output, fileSave)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 23:32:58
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.052413
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 3, for S/N=7.015
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 30.0
Flux in a centered 1xFWHM circular aperture = 1.948
Central pixel S/N = 8.027
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 6.612
Max S/N (shifting the aperture center) = 9.537
stddev S/N (shifting the aperture center) = 1.650
../_images/tutorials_07_ifs_psfsub_fm_planets_215_1.png
../_images/tutorials_07_ifs_psfsub_fm_planets_215_2.png 
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 23:33:00
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [83.75877533353938,29.974810067878558], running ...
Step | flux    | chi2r
1/30   0.100   0.081
2/30   0.149   0.080
3/30   0.221   0.079
4/30   0.329   0.077
5/30   0.489   0.075
6/30   0.728   0.071
7/30   1.083   0.066
8/30   1.610   0.059
9/30   2.395   0.050
10/30   3.562   0.040
11/30   5.298   0.032
12/30   7.880   0.035
13/30   11.721   0.072
14/30   17.433   0.211
../_images/tutorials_07_ifs_psfsub_fm_planets_215_4.png 
Planet 0: preliminary position guess: (r, theta)=(38.4, 323.2)
Planet 0: preliminary flux guess: 5.3
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 24, nfev: 55, chi2r: 0.030479820234548226
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f)=(38.422, 323.182, 5.986) at
          (X,Y)=(83.76, 29.97)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:05.709994
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-04 23:33:05
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
        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.00 and 0.10 respectively.
Beginning emcee Ensemble sampler...
emcee Ensemble sampler successful

Start of the MCMC run ...
Step  |  Duration/step (sec)  |  Remaining Estimated Time (sec)
0               3.47809                 1735.56816
1               3.54162                 1763.72676
2               3.56924                 1773.91129
3               3.57911                 1775.24054
4               3.53748                 1751.05359
5               3.56513                 1761.17521
6               3.59093                 1770.32800
7               3.59911                 1770.76015
8               3.56543                 1750.62564
9               3.56631                 1747.49141
10              3.57485                 1748.10263
11              3.39883                 1658.62953
12              3.50620                 1707.52086
13              3.55256                 1726.54659
14              3.43421                 1665.59331
15              3.35669                 1624.63699
16              3.35914                 1622.46703
17              3.25477                 1568.79914
18              3.44416                 1656.64336
19              3.52268                 1690.88400
20              3.35348                 1606.31596
21              3.39113                 1620.95871
22              3.41184                 1627.44959
23              3.21538                 1530.51898
24              3.24365                 1540.73375
25              3.37631                 1600.36999
26              3.47933                 1645.72167
27              3.28789                 1551.88172
28              3.30518                 1556.73743
29              3.49941                 1644.72223
30              3.31703                 1555.68613
31              3.40476                 1593.42815
32              3.50548                 1637.06103
33              3.40064                 1584.69964
34              3.36067                 1562.71155
35              3.52626                 1636.18371
36              3.42797                 1587.15057
37              3.38784                 1565.18347
38              3.38397                 1560.01109
39              3.41991                 1573.15814
40              3.52405                 1617.53757
41              3.26755                 1496.53561
42              3.38184                 1545.49997
43              3.41352                 1556.56512
44              3.20909                 1460.13777
45              3.45107                 1566.78760
46              3.40573                 1542.79614
47              3.10866                 1405.11342
48              3.27754                 1478.17189
49              3.15829                 1421.22870
50              3.25500                 1461.49455

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_6.png 
51              3.37048                 1509.97594
52              3.24713                 1451.46622
53              3.13319                 1397.40051
54              3.29470                 1466.14194
55              3.21198                 1426.12045
56              3.19556                 1415.63087
57              3.18334                 1407.03540
58              2.84094                 1252.85322
59              3.14935                 1385.71576
60              3.07052                 1347.95872
61              3.24476                 1421.20444
62              3.26626                 1427.35431
63              3.23666                 1411.18376
64              2.94819                 1282.46309
65              3.10262                 1346.53665
66              3.28214                 1421.16705
67              3.21433                 1388.58970
68              3.34086                 1439.90937
69              3.22403                 1386.33204
70              3.30994                 1419.96426
71              3.15832                 1351.76182
72              3.23329                 1380.61355
73              3.21582                 1369.93932
74              3.20695                 1362.95418
75              3.20002                 1356.80763
76              3.21823                 1361.30960
77              3.35177                 1414.44694
78              3.26758                 1375.64950
79              3.25456                 1366.91352
80              3.39811                 1423.80851
81              3.23085                 1350.49655
82              3.31778                 1383.51509
83              3.37380                 1403.50288
84              3.26490                 1354.93225
85              3.26864                 1353.21489
86              3.24233                 1339.08312
87              3.15044                 1297.98046
88              3.05068                 1253.83071
89              3.27804                 1343.99681
90              3.34270                 1367.16307
91              3.25281                 1327.14689
92              3.64683                 1484.25818
93              3.40264                 1381.47143
94              3.38121                 1369.39086
95              3.30169                 1333.88195
96              3.29015                 1325.93085
97              3.25212                 1307.35103
98              3.11156                 1247.73516
99              3.18515                 1274.05960
100             3.25925                 1300.44075

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_8.png 
101             3.16939                 1261.41523
102             3.35320                 1331.21961
103             3.04316                 1205.08938
104             3.30938                 1307.20668
105             3.22802                 1271.84146
106             3.29518                 1295.00692
107             3.27320                 1283.09362
108             3.26464                 1276.47424
109             3.22700                 1258.53039
110             3.16030                 1229.35826
111             2.94485                 1142.60296
112             3.16645                 1225.41692
113             2.97331                 1147.69805
114             3.23068                 1243.81141
115             3.16268                 1214.46989
116             3.14318                 1203.83947
117             3.22305                 1231.20472
118             3.34750                 1275.39598
119             3.20051                 1216.19380
120             3.15401                 1195.36827
121             3.12741                 1182.15947
122             3.19434                 1204.26543
123             3.34615                 1258.15315
124             3.27010                 1226.28900
125             3.28646                 1229.13679
126             3.16084                 1178.99183
127             3.13848                 1167.51344
128             3.21954                 1194.44934
129             3.20780                 1186.88674
130             3.15889                 1165.63115
131             3.26631                 1202.00282
132             3.21896                 1181.35905
133             3.11868                 1141.43542
134             3.26778                 1192.74043
135             3.26742                 1189.34088
136             3.22504                 1170.68771
137             3.15220                 1141.09821
138             3.19289                 1152.63293
139             3.43362                 1236.10428
140             3.14692                 1129.74320
141             3.17623                 1137.09177
142             3.30993                 1181.64394
143             3.06021                 1089.43654
144             3.15303                 1119.32423
145             3.27154                 1158.12445
146             3.35046                 1182.71273
147             3.38557                 1191.71958
148             3.25451                 1142.33161
149             3.10742                 1087.59700
150             3.21516                 1122.09119

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_10.png 
151             3.22027                 1120.65361
152             3.10274                 1076.65147
153             3.24997                 1124.49100
154             3.25184                 1121.88549
155             3.24569                 1116.51839
156             3.18141                 1091.22500
157             3.16482                 1082.36776
158             3.21492                 1096.28874
159             3.25624                 1107.12262
160             3.22056                 1091.76848
161             3.15461                 1066.25784
162             3.11503                 1049.76376
163             3.23575                 1087.21267
164             3.23549                 1083.88949
165             3.07234                 1026.16323
166             3.26970                 1088.81010
167             3.12652                 1038.00331
168             3.31556                 1097.45003
169             2.95299                 974.48571
170             3.17184                 1043.53668
171             3.03849                 996.62505
172             3.18132                 1040.29229
173             3.29925                 1075.55680
174             3.26096                 1059.81298
175             3.03235                 982.48172
176             3.26329                 1054.04364
177             3.15094                 1014.60397
178             3.15680                 1013.33152
179             3.30095                 1056.30464
180             3.41898                 1090.65526
181             3.26507                 1038.29226
182             3.31988                 1052.40228
183             3.26105                 1030.49148
184             3.40993                 1074.12732
185             3.41030                 1070.83546
186             3.11075                 973.66350
187             3.17277                 989.90362
188             3.29908                 1026.01450
189             3.29292                 1020.80489
190             3.27132                 1010.83912
191             3.05238                 940.13150
192             3.11632                 956.70963
193             3.08949                 945.38486
194             3.22589                 983.89554
195             3.37860                 1027.09562
196             3.32559                 1007.65347
197             3.28625                 992.44780
198             3.32438                 1000.63988
199             2.96768                 890.30550
200             3.20728                 958.97642

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_12.png 
Auto-corr tau/N = [0.05552478]
tau/N <= 0.02 = [False]

201             3.30833                 985.88174
202             3.18899                 947.12973
203             3.06807                 908.15020
204             3.17291                 936.00904
205             3.10636                 913.26866
206             3.01203                 882.52508
207             3.19787                 933.77892
208             3.34239                 972.63665
209             3.36180                 974.92200
210             3.29214                 951.42701
211             3.33859                 961.51306
212             3.22409                 925.31498
213             3.29306                 941.81402
214             3.16160                 901.05685
215             3.11457                 884.53731
216             3.09325                 875.39088
217             3.35819                 947.00986
218             3.11049                 874.04713
219             3.26494                 914.18432
220             3.27158                 912.77138
221             3.10809                 864.04985
222             3.20433                 887.59858
223             3.22140                 889.10750
224             3.05652                 840.54327
225             3.31937                 909.50711
226             3.37083                 920.23741
227             3.35700                 913.10318
228             3.16520                 857.77001
229             3.26564                 881.72226
230             3.21057                 863.64387
231             3.35398                 898.86664
232             3.28969                 878.34643
233             3.20135                 851.56016
234             3.25981                 863.84912
235             3.29273                 869.28204
236             3.11830                 820.11369
237             3.20796                 840.48604
238             3.36063                 877.12365
239             3.06559                 797.05366
240             3.32614                 861.47155
241             3.17044                 817.97300
242             3.12649                 803.50921
243             3.22185                 824.79283
244             3.02123                 770.41289
245             3.03707                 771.41451
246             3.20734                 811.45803
247             3.15675                 795.50024
248             3.15980                 793.10955
249             3.26587                 816.46825
250             3.18195                 792.30630

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_14.png 
Auto-corr tau/N = [0.05248042]
tau/N <= 0.02 = [False]

251             3.60909                 895.05407
252             3.18838                 787.52961
253             3.05312                 751.06850
254             3.25482                 797.43188
255             3.15482                 769.77510
256             3.23075                 785.07249
257             3.24331                 784.88199
258             3.21443                 774.67787
259             3.17194                 761.26632
260             3.31989                 793.45275
261             3.30177                 785.82031
262             3.01880                 715.45536
263             3.20989                 757.53498
264             3.16967                 744.87315
265             3.25547                 761.78021
266             3.16208                 736.76464
267             3.34522                 776.09104
268             3.29601                 761.37900
269             3.34110                 768.45254
270             3.15868                 723.33749
271             3.21232                 732.40850
272             3.43062                 778.74983
273             3.34666                 756.34561
274             3.33864                 751.19445
275             3.21529                 720.22563
276             3.24921                 724.57383
277             3.26218                 724.20374
278             3.26935                 722.52613
279             3.24201                 713.24308
280             3.38235                 740.73553
281             3.26464                 711.69261
282             3.31855                 720.12448
283             3.29479                 711.67550
284             3.32951                 715.84401
285             3.33020                 712.66280
286             3.22962                 687.90991
287             3.37880                 716.30539
288             3.35429                 707.75456
289             3.28975                 690.84750
290             3.17855                 664.31653
291             3.14645                 654.46181
292             3.07328                 636.16917
293             3.26768                 673.14126
294             3.14200                 644.11021
295             3.24864                 662.72297
296             3.20666                 650.95096
297             3.11193                 628.61026
298             3.24641                 652.52741
299             3.15151                 630.30260
300             3.19737                 636.27603

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_16.png 
Auto-corr tau/N = [0.05102172]
tau/N <= 0.02 = [False]

301             3.28032                 649.50395
302             3.25318                 640.87685
303             3.33528                 653.71547
304             3.23511                 630.84684
305             3.08401                 598.29736
306             3.07061                 592.62715
307             3.06763                 588.98534
308             3.07704                 587.71521
309             3.16415                 601.18812
310             3.30228                 624.13111
311             3.21719                 604.83191
312             3.05661                 571.58607
313             3.23822                 602.30948
314             3.30636                 611.67678
315             3.23007                 594.33288
316             3.19307                 584.33126
317             3.28137                 597.20843
318             3.30480                 598.16826
319             3.24846                 584.72352
320             3.27619                 586.43747
321             3.30036                 587.46479
322             3.26374                 577.68233
323             3.28715                 578.53787
324             3.09854                 542.24468
325             3.28268                 571.18684
326             3.30488                 571.74476
327             3.24792                 558.64241
328             3.25005                 555.75838
329             3.21187                 546.01739
330             3.23423                 546.58470
331             2.94072                 494.04096
332             3.16631                 528.77410
333             3.28771                 545.75969
334             3.17369                 523.65852
335             3.27966                 537.86457
336             3.20092                 521.74980
337             3.17681                 514.64354
338             3.43955                 553.76771
339             3.25906                 521.44896
340             3.18736                 506.78976
341             3.20673                 506.66350
342             3.41306                 535.85011
343             3.10903                 485.00899
344             3.14652                 487.71122
345             3.06829                 472.51697
346             3.08303                 471.70359
347             3.23211                 491.28057
348             3.12292                 471.56167
349             2.99031                 448.54665
350             3.23232                 481.61523

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_18.png 
Auto-corr tau/N = [0.0473208]
tau/N <= 0.02 = [False]

351             3.20762                 474.72717
352             2.79757                 411.24220
353             3.11444                 454.70824
354             3.26143                 472.90691
355             3.10719                 447.43594
356             3.04172                 434.96610
357             3.14800                 447.01543
358             3.21053                 452.68417
359             3.19389                 447.14446
360             3.23214                 449.26732
361             3.13200                 432.21669
362             3.21975                 441.10589
363             3.52880                 479.91680
364             3.07911                 415.68012
365             3.24319                 434.58746
366             3.23309                 430.00110
367             3.20237                 422.71231
368             3.04260                 398.58099
369             3.29259                 428.03709
370             3.27764                 422.81556
371             3.23461                 414.02970
372             3.32605                 422.40771
373             3.12142                 393.29854
374             3.10860                 388.57438
375             3.20931                 397.95407
376             3.16374                 389.13990
377             3.24977                 396.47170
378             3.34674                 404.95578
379             3.28110                 393.73164
380             3.09922                 368.80766
381             3.09285                 364.95618
382             3.15241                 368.83232
383             3.13178                 363.28590
384             3.40766                 391.88136
385             3.19583                 364.32428
386             3.23556                 365.61783
387             3.35435                 375.68731
388             3.48659                 387.01138
389             3.41819                 376.00134
390             3.28911                 358.51353
391             3.44992                 372.59190
392             3.32942                 356.24826
393             3.40583                 361.01798
394             3.40019                 357.02006
395             3.19062                 331.82417
396             3.14641                 324.08044
397             3.31008                 337.62816
398             3.22436                 325.66006
399             3.60796                 360.79570
400             3.17251                 314.07809

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_20.png 
Auto-corr tau/N = [0.04277931]
tau/N <= 0.02 = [False]

401             3.17008                 310.66774
402             3.21624                 311.97499
403             3.25553                 312.53069
404             3.18553                 302.62554
405             3.14325                 295.46531
406             3.07774                 286.23001
407             3.06590                 282.06289
408             3.26145                 296.79186
409             3.06308                 275.67729
410             3.30491                 294.13681
411             2.94921                 259.53013
412             3.12289                 271.69126
413             3.06270                 263.39229
414             3.20000                 272.00000
415             3.38106                 284.00870
416             3.27517                 271.83903
417             3.26391                 267.64054
418             3.27641                 265.38945
419             3.16116                 252.89256
420             3.13655                 247.78761
421             3.18884                 248.72936
422             3.05119                 234.94171
423             3.13500                 238.26000
424             3.06861                 230.14575
425             3.20547                 237.20456
426             3.01137                 219.82994
427             3.13004                 225.36324
428             3.13977                 222.92367
429             3.24261                 226.98235
430             3.13129                 216.05901
431             3.00688                 204.46764
432             3.18640                 213.48860
433             3.18154                 209.98138
434             3.13353                 203.67932
435             3.06537                 196.18381
436             3.19325                 201.17500
437             3.28757                 203.82909
438             3.10778                 189.57482
439             3.19320                 191.59182
440             3.22055                 190.01269
441             3.22309                 186.93922
442             3.19810                 182.29193
443             3.21578                 180.08379
444             3.29389                 181.16378
445             3.32446                 179.52068
446             3.11152                 164.91067
447             3.15674                 164.15074
448             3.14663                 160.47839
449             3.05998                 152.99900
450             3.24347                 158.93013

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_22.png 
Auto-corr tau/N = [0.03904012]
tau/N <= 0.02 = [False]

451             3.49893                 167.94878
452             3.11547                 146.42695
453             3.23002                 148.58087
454             3.02799                 136.25950
455             3.21790                 141.58764
456             3.24703                 139.62246
457             3.14977                 132.29042
458             3.14687                 129.02167
459             3.25036                 130.01456
460             3.22131                 125.63105
461             3.04470                 115.69856
462             3.19686                 118.28382
463             3.34044                 120.25591
464             3.13578                 109.75234
465             3.28993                 111.85769
466             3.24694                 107.14892
467             3.18162                 101.81181
468             3.19404                 99.01512
469             3.06844                 92.05323
470             2.83036                 82.08053
471             3.12161                 87.40508
472             3.17282                 85.66619
473             2.82664                 73.49274
474             3.24429                 81.10730
475             3.28950                 78.94788
476             2.99987                 68.99694
477             3.02706                 66.59532
478             2.96212                 62.20456
479             3.08509                 61.70184
480             3.18573                 60.52893
481             3.06821                 55.22787
482             3.06782                 52.15294
483             3.12170                 49.94722
484             3.12633                 46.89487
485             3.04772                 42.66815
486             3.00621                 39.08074
487             3.13239                 37.58863
488             3.16605                 34.82657
489             3.24451                 32.44512
490             3.24593                 29.21340
491             3.07867                 24.62935
492             3.34766                 23.43365
493             3.05389                 18.32333
494             3.18970                 15.94851
495             3.12846                 12.51385
496             3.24037                 9.72112
497             3.14605                 6.29211
498             3.21718                 3.21718
499             3.06088                 0.00000
We have reached the limit # of steps without convergence
Running time:  0:26:59.034150
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-05 00:00:04
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.080102
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 5, for S/N=16.898
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 30.0
Flux in a centered 1xFWHM circular aperture = 3.720
Central pixel S/N = 13.660
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 10.367
Max S/N (shifting the aperture center) = 15.549
stddev S/N (shifting the aperture center) = 3.361
../_images/tutorials_07_ifs_psfsub_fm_planets_215_24.png
../_images/tutorials_07_ifs_psfsub_fm_planets_215_25.png 
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-05 00:00:06
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [83.75877533353938,29.974810067878558], running ...
Step | flux    | chi2r
1/30   0.100   0.205
2/30   0.149   0.204
3/30   0.221   0.203
4/30   0.329   0.202
5/30   0.489   0.199
6/30   0.728   0.196
7/30   1.083   0.190
8/30   1.610   0.182
9/30   2.395   0.169
10/30   3.562   0.153
11/30   5.298   0.130
12/30   7.880   0.100
13/30   11.721   0.064
14/30   17.433   0.031
15/30   25.929   0.042
16/30   38.566   0.176
17/30   57.362   0.402
../_images/tutorials_07_ifs_psfsub_fm_planets_215_27.png 
Planet 0: preliminary position guess: (r, theta)=(38.4, 323.2)
Planet 0: preliminary flux guess: 17.4
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 27, nfev: 59, chi2r: 0.02695430941500906
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f)=(38.422, 323.182, 20.304) at
          (X,Y)=(83.76, 29.97)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:06.400024
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-05 00:00:12
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
        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.00 and 0.11 respectively.
Beginning emcee Ensemble sampler...
emcee Ensemble sampler successful

Start of the MCMC run ...
Step  |  Duration/step (sec)  |  Remaining Estimated Time (sec)
0               3.39102                 1692.11698
1               3.46186                 1724.00578
2               3.47339                 1726.27483
3               3.49464                 1733.34392
4               3.52650                 1745.61503
5               3.45711                 1707.81036
6               3.45780                 1704.69589
7               3.45956                 1702.10155
8               3.46334                 1700.50239
9               3.65807                 1792.45626
10              3.45305                 1688.54292
11              3.49467                 1705.40091
12              3.49205                 1700.62884
13              3.45106                 1677.21467
14              3.50978                 1702.24378
15              3.48760                 1688.00082
16              3.46787                 1674.98314
17              3.49829                 1686.17771
18              3.49024                 1678.80304
19              3.49207                 1676.19552
20              3.49269                 1673.00043
21              3.50853                 1677.07686
22              3.48787                 1663.71590
23              3.49351                 1662.90933
24              3.48901                 1657.28165
25              3.45793                 1639.05787
26              3.49283                 1652.11001
27              3.47679                 1641.04299
28              3.47289                 1635.73025
29              3.58137                 1683.24531
30              3.54799                 1664.00825
31              3.52748                 1650.86298
32              3.52740                 1647.29440
33              3.48890                 1625.82787
34              3.49818                 1626.65231
35              3.45861                 1604.79643
36              3.46801                 1605.69094
37              3.45883                 1597.97900
38              3.45224                 1591.48126
39              3.46269                 1592.83602
40              3.44785                 1582.56177
41              3.45428                 1582.05932
42              3.45338                 1578.19649
43              3.58144                 1633.13892
44              3.46022                 1574.40101
45              3.45720                 1569.56698
46              3.51807                 1593.68707
47              3.55046                 1604.80656
48              3.49695                 1577.12220
49              3.47831                 1565.23905
50              3.47501                 1560.27949

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_29.png 
51              3.62630                 1624.58374
52              3.50838                 1568.24363
53              3.45686                 1541.76090
54              3.45681                 1538.28134
55              3.46203                 1537.13999
56              3.46698                 1535.87170
57              3.46032                 1529.46100
58              3.46636                 1528.66388
59              3.46243                 1523.46744
60              3.47415                 1525.15361
61              3.75600                 1645.12712
62              3.47404                 1518.15373
63              3.49081                 1521.99490
64              3.49850                 1521.84750
65              3.41319                 1481.32316
66              3.35843                 1454.19976
67              3.42634                 1480.17715
68              3.45143                 1487.56805
69              3.38841                 1457.01630
70              3.43694                 1474.44683
71              3.45434                 1478.45624
72              3.40707                 1454.82060
73              3.42422                 1458.71942
74              3.35964                 1427.84912
75              3.45453                 1464.71987
76              3.46903                 1467.39842
77              3.43067                 1447.74485
78              3.58528                 1509.40078
79              3.47526                 1459.60836
80              3.43301                 1438.43245
81              3.42732                 1432.62101
82              3.38198                 1410.28358
83              3.36272                 1398.89152
84              3.36370                 1395.93384
85              3.36242                 1392.04105
86              3.37464                 1393.72591
87              3.36309                 1385.59184
88              3.36479                 1382.92992
89              3.44636                 1413.00719
90              3.44825                 1410.33507
91              3.43022                 1399.52772
92              3.41404                 1389.51591
93              3.36628                 1366.70887
94              3.42343                 1386.48955
95              3.46398                 1399.44954
96              3.38846                 1365.54817
97              3.38337                 1360.11313
98              3.42991                 1375.39311
99              3.37743                 1350.97200
100             3.37263                 1345.68136

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_31.png 
101             3.50522                 1395.07915
102             3.36391                 1335.47148
103             3.36861                 1333.96798
104             3.37058                 1331.37949
105             3.37837                 1331.07699
106             3.37285                 1325.53005
107             3.36542                 1319.24621
108             3.38019                 1321.65351
109             3.37294                 1315.44777
110             3.37292                 1312.06744
111             3.38619                 1313.84056
112             3.37971                 1307.94622
113             3.49099                 1347.52137
114             3.38830                 1304.49435
115             3.38411                 1299.49709
116             3.40409                 1303.76838
117             3.37019                 1287.41182
118             3.36841                 1283.36535
119             3.38335                 1285.67262
120             3.37919                 1280.71149
121             3.37827                 1276.98455
122             3.38474                 1276.04736
123             3.37169                 1267.75394
124             3.37530                 1265.73600
125             3.37413                 1261.92387
126             3.36866                 1256.51167
127             3.36595                 1252.13489
128             3.36509                 1248.44691
129             3.37261                 1247.86644
130             3.37905                 1246.87093
131             3.37987                 1243.79142
132             3.37561                 1238.84777
133             3.36782                 1232.62102
134             3.40473                 1242.72609
135             3.37198                 1227.40145
136             3.37961                 1226.79952
137             3.36634                 1218.61653
138             3.36892                 1216.18048
139             3.37392                 1214.60976
140             3.36971                 1209.72553
141             3.37006                 1206.48220
142             3.37411                 1204.55549
143             3.37690                 1202.17462
144             3.36517                 1194.63535
145             3.36668                 1191.80437
146             3.37212                 1190.35836
147             3.37641                 1188.49456
148             3.40758                 1196.05918
149             3.38808                 1185.82625
150             3.36704                 1175.09591

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_33.png 
151             3.52671                 1227.29473
152             3.36980                 1169.32025
153             3.36237                 1163.37829
154             3.36511                 1160.96260
155             3.36028                 1155.93735
156             3.36209                 1153.19584
157             3.37359                 1153.76846
158             3.37824                 1151.97848
159             3.36034                 1142.51628
160             3.36438                 1140.52448
161             3.37564                 1140.96700
162             3.37358                 1136.89646
163             3.37620                 1134.40186
164             3.37794                 1131.61023
165             3.37435                 1127.03390
166             3.37063                 1122.42012
167             3.37252                 1119.67564
168             3.36939                 1115.26776
169             3.38450                 1116.88599
170             3.37952                 1111.86307
171             3.36069                 1102.30763
172             3.36249                 1099.53554
173             3.36954                 1098.47134
174             3.36312                 1093.01465
175             3.36451                 1090.10124
176             3.37008                 1088.53455
177             3.37686                 1087.34795
178             3.37315                 1082.78019
179             3.36897                 1078.07136
180             3.36878                 1074.64050
181             3.37770                 1074.10987
182             3.37085                 1068.55977
183             3.37268                 1065.76720
184             3.44911                 1086.47091
185             3.37308                 1059.14712
186             3.38339                 1059.00013
187             3.38720                 1056.80671
188             3.36988                 1048.03206
189             3.36854                 1044.24616
190             3.40848                 1053.22032
191             3.36567                 1036.62759
192             3.36734                 1033.77430
193             3.37774                 1033.58844
194             3.36467                 1026.22374
195             3.36728                 1023.65434
196             3.37673                 1023.14889
197             3.37449                 1019.09538
198             3.36357                 1012.43577
199             3.36153                 1008.45900
200             3.37769                 1009.92961

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_35.png 
Auto-corr tau/N = [0.0628609]
tau/N <= 0.02 = [False]

201             3.52598                 1050.74323
202             3.37882                 1003.50835
203             3.36756                 996.79894
204             3.37247                 994.87865
205             3.42193                 1006.04860
206             3.37695                 989.44518
207             3.37934                 986.76728
208             3.37544                 982.25275
209             3.36451                 975.70732
210             3.37316                 974.84411
211             3.36896                 970.25962
212             3.37390                 968.30786
213             3.38426                 967.89693
214             3.37922                 963.07856
215             3.37125                 957.43472
216             3.37683                 955.64232
217             3.36732                 949.58396
218             3.36768                 946.31724
219             3.38181                 946.90764
220             3.41065                 951.57079
221             3.38270                 940.39060
222             3.40496                 943.17309
223             3.38128                 933.23411
224             3.37825                 929.01875
225             3.37956                 926.00054
226             3.39586                 927.07087
227             3.37356                 917.60778
228             3.36829                 912.80767
229             3.36728                 909.16479
230             3.37088                 906.76591
231             3.36649                 902.22012
232             3.39615                 906.77285
233             3.39196                 902.26189
234             3.36703                 892.26216
235             3.36593                 888.60578
236             3.37709                 888.17362
237             3.37161                 883.36261
238             3.36348                 877.86932
239             3.39631                 883.04164
240             3.46897                 898.46297
241             3.38901                 874.36355
242             3.36534                 864.89264
243             3.39629                 869.45126
244             3.41732                 871.41609
245             3.38007                 858.53880
246             3.37158                 853.00949
247             3.37369                 850.17089
248             3.36817                 845.40992
249             3.36542                 841.35400
250             3.36527                 837.95323

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_37.png 
Auto-corr tau/N = [0.05818584]
tau/N <= 0.02 = [False]

251             3.53661                 877.08052
252             3.41670                 843.92490
253             3.36992                 829.00007
254             3.36927                 825.47189
255             3.44945                 841.66556
256             3.41234                 829.19813
257             3.41560                 826.57593
258             3.45221                 831.98237
259             3.36591                 807.81744
260             3.37182                 805.86617
261             3.37970                 804.36931
262             3.38019                 801.10550
263             3.37462                 796.40938
264             3.60018                 846.04254
265             3.45231                 807.84031
266             3.39948                 792.07791
267             3.41562                 792.42500
268             3.39255                 783.67882
269             3.40104                 782.23966
270             3.53186                 808.79548
271             3.51948                 802.44144
272             3.48115                 790.21992
273             3.40835                 770.28665
274             3.42555                 770.74785
275             3.41957                 765.98435
276             3.40215                 758.68012
277             3.36880                 747.87404
278             3.37641                 746.18727
279             3.36559                 740.43002
280             3.37122                 738.29762
281             3.36565                 733.71192
282             3.38693                 734.96359
283             3.41304                 737.21772
284             3.40484                 732.04125
285             3.37685                 722.64526
286             3.39141                 722.37033
287             3.38287                 717.16802
288             3.39831                 717.04299
289             3.42192                 718.60404
290             3.45042                 721.13778
291             3.41464                 710.24470
292             3.38867                 701.45490
293             3.46999                 714.81753
294             3.42367                 701.85317
295             3.39743                 693.07572
296             3.39513                 689.21037
297             3.40663                 688.13886
298             3.37158                 677.68838
299             3.37344                 674.68820
300             3.36831                 670.29309

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_39.png 
Auto-corr tau/N = [0.05007516]
tau/N <= 0.02 = [False]

301             3.57329                 707.51083
302             3.54739                 698.83544
303             3.38305                 663.07858
304             3.37418                 657.96549
305             3.38810                 657.29198
306             3.42966                 661.92419
307             3.41394                 655.47686
308             3.36583                 642.87410
309             3.38394                 642.94879
310             3.37883                 638.59868
311             3.53982                 665.48654
312             3.36521                 629.29502
313             3.39152                 630.82198
314             3.37164                 623.75396
315             3.37125                 620.30908
316             3.37050                 616.80113
317             3.39037                 617.04679
318             3.45578                 625.49582
319             3.40952                 613.71342
320             3.39133                 607.04789
321             3.36490                 598.95238
322             3.36375                 595.38393
323             3.38136                 595.11883
324             3.36727                 589.27278
325             3.44534                 599.48986
326             3.38205                 585.09551
327             3.36170                 578.21274
328             3.39617                 580.74490
329             3.36733                 572.44644
330             3.35771                 567.45333
331             3.36456                 565.24608
332             3.35973                 561.07474
333             3.36333                 558.31228
334             3.36634                 555.44561
335             3.36574                 551.98103
336             3.35847                 547.43110
337             3.36132                 544.53352
338             3.36197                 541.27653
339             3.36437                 538.29936
340             3.35791                 533.90705
341             3.36449                 531.58879
342             3.35803                 527.21071
343             3.32134                 518.12888
344             3.40901                 528.39577
345             3.39001                 522.06092
346             3.43809                 526.02808
347             3.40078                 516.91902
348             3.37218                 509.19948
349             3.36024                 504.03615
350             3.35958                 500.57712

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_41.png 
Auto-corr tau/N = [0.0468243]
tau/N <= 0.02 = [False]

351             3.72194                 550.84771
352             3.42347                 503.25024
353             3.38095                 493.61826
354             3.37273                 489.04556
355             3.36955                 485.21578
356             3.37449                 482.55150
357             3.37212                 478.84076
358             3.37827                 476.33635
359             3.36330                 470.86256
360             3.37225                 468.74233
361             3.45326                 476.54947
362             3.36916                 461.57451
363             3.36965                 458.27267
364             3.38091                 456.42258
365             3.35996                 450.23451
366             3.38602                 450.34106
367             3.46579                 457.48402
368             3.41362                 447.18435
369             3.37946                 439.32928
370             3.39695                 438.20616
371             3.43837                 440.11149
372             3.39206                 430.79137
373             3.40107                 428.53520
374             3.38476                 423.09438
375             3.38093                 419.23544
376             3.39069                 417.05450
377             3.37684                 411.97460
378             3.38267                 409.30331
379             3.36583                 403.89936
380             3.35803                 399.60569
381             3.36447                 397.00781
382             3.38698                 396.27666
383             3.37425                 391.41242
384             3.36418                 386.88093
385             3.35950                 382.98357
386             3.36158                 379.85865
387             3.40867                 381.77149
388             3.37623                 374.76142
389             3.36253                 369.87863
390             3.49863                 381.35100
391             3.36916                 363.86939
392             3.46552                 370.81096
393             3.45364                 366.08542
394             3.48620                 366.05100
395             3.38319                 351.85145
396             3.43580                 353.88719
397             3.45895                 352.81321
398             3.42730                 346.15710
399             3.43497                 343.49730
400             3.44689                 341.24171

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_43.png 
Auto-corr tau/N = [0.04295486]
tau/N <= 0.02 = [False]

401             3.80881                 373.26348
402             3.36186                 326.10071
403             3.36931                 323.45357
404             3.37133                 320.27663
405             3.36821                 316.61136
406             3.44276                 320.17705
407             3.37512                 310.51122
408             3.42603                 311.76864
409             3.37240                 303.51564
410             3.36343                 299.34509
411             3.40746                 299.85674
412             3.37679                 293.78038
413             3.40218                 292.58791
414             3.36168                 285.74280
415             3.38924                 284.69624
416             3.42396                 284.18868
417             3.45509                 283.31713
418             3.39556                 275.04012
419             3.38693                 270.95448
420             3.43017                 270.98327
421             3.44792                 268.93737
422             3.51367                 270.55290
423             3.39060                 257.68598
424             3.39589                 254.69205
425             3.41848                 252.96767
426             3.44704                 251.63363
427             3.40573                 245.21292
428             3.36802                 239.12977
429             3.42246                 239.57192
430             3.37448                 232.83926
431             3.37252                 229.33129
432             3.43152                 229.91151
433             3.37985                 223.06984
434             3.44186                 223.72123
435             3.48939                 223.32083
436             3.45748                 217.82099
437             3.41743                 211.88078
438             3.40855                 207.92167
439             3.39074                 203.44464
440             3.39503                 200.30689
441             3.43106                 199.00154
442             3.40154                 193.88767
443             3.37384                 188.93487
444             3.37134                 185.42392
445             3.41548                 184.43587
446             3.39092                 179.71892
447             3.42593                 178.14862
448             3.37941                 172.35006
449             3.42984                 171.49190
450             3.42688                 167.91707

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_45.png 
Auto-corr tau/N = [0.03980261]
tau/N <= 0.02 = [False]

451             3.52138                 169.02638
452             3.45667                 162.46354
453             3.36849                 154.95077
454             3.36058                 151.22628
455             3.45686                 152.10184
456             3.39073                 145.80152
457             3.42521                 143.85865
458             3.38245                 138.68057
459             3.36840                 134.73584
460             3.36430                 131.20790
461             3.43270                 130.44245
462             3.42481                 126.71801
463             3.38861                 121.98978
464             3.38178                 118.36216
465             3.38719                 115.16449
466             3.37644                 111.42252
467             3.37149                 107.88758
468             3.39642                 105.28905
469             3.40392                 102.11763
470             3.39758                 98.52988
471             3.40207                 95.25790
472             3.38968                 91.52125
473             3.44038                 89.44978
474             3.42506                 85.62655
475             3.46207                 83.08963
476             3.43163                 78.92751
477             3.41604                 75.15299
478             3.40643                 71.53507
479             3.36077                 67.21546
480             3.36053                 63.85009
481             3.37116                 60.68086
482             3.35920                 57.10640
483             3.35798                 53.72773
484             3.40641                 51.09616
485             3.38474                 47.38633
486             3.36130                 43.69693
487             3.40021                 40.80256
488             3.37342                 37.10764
489             3.40441                 34.04412
490             3.45881                 31.12930
491             3.44787                 27.58292
492             3.40154                 23.81077
493             3.46087                 20.76521
494             3.47075                 17.35374
495             3.44011                 13.76043
496             3.43981                 10.31944
497             3.43660                 6.87319
498             3.46035                 3.46035
499             3.40095                 0.00000
We have reached the limit # of steps without convergence
Running time:  0:28:28.713779
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-05 00:28:41
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Done SVD/PCA with numpy SVD (LAPACK)
Running time:  0:00:00.024079
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Number of steps 11
Optimal number of PCs = 2, for S/N=23.064
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Coords of chosen px (X,Y) = 83.8, 30.0
Flux in a centered 1xFWHM circular aperture = 7.519
Central pixel S/N = 17.995
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Inside a centered 1xFWHM circular aperture:
Mean S/N (shifting the aperture center) = 14.409
Max S/N (shifting the aperture center) = 19.392
stddev S/N (shifting the aperture center) = 3.291
../_images/tutorials_07_ifs_psfsub_fm_planets_215_47.png
../_images/tutorials_07_ifs_psfsub_fm_planets_215_48.png 
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-05 00:28:42
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
             Planet 0
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Planet 0: flux estimation at the position [83.75877533353938,29.974810067878558], running ...
Step | flux    | chi2r
1/30   0.100   0.627
2/30   0.149   0.625
3/30   0.221   0.622
4/30   0.329   0.618
5/30   0.489   0.612
6/30   0.728   0.602
7/30   1.083   0.588
8/30   1.610   0.567
9/30   2.395   0.537
10/30   3.562   0.493
11/30   5.298   0.431
12/30   7.880   0.349
13/30   11.721   0.229
14/30   17.433   0.107
15/30   25.929   0.032
16/30   38.566   0.201
17/30   57.362   0.902
18/30   85.317   2.616
../_images/tutorials_07_ifs_psfsub_fm_planets_215_50.png 
Planet 0: preliminary position guess: (r, theta)=(38.4, 323.2)
Planet 0: preliminary flux guess: 25.9
Planet 0: Simplex Nelder-Mead minimization, running ...
Planet 0: Success: True, nit: 22, nfev: 49, chi2r: 0.031329481433536005
message: Optimization terminated successfully.
Planet 0: simplex result: (r, theta, f)=(38.422, 323.182, 26.056) at
          (X,Y)=(83.76, 29.97)

 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
DONE !
 ――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Running time:  0:00:05.363499
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Starting time: 2022-05-05 00:28:48
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
        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.00 and 0.10 respectively.
Beginning emcee Ensemble sampler...
emcee Ensemble sampler successful

Start of the MCMC run ...
Step  |  Duration/step (sec)  |  Remaining Estimated Time (sec)
0               3.45894                 1726.01281
1               3.44648                 1716.34654
2               3.39917                 1689.38600
3               3.35724                 1665.18955
4               3.35719                 1661.80856
5               3.42947                 1694.15620
6               3.39981                 1676.10534
7               3.36700                 1656.56498
8               3.35496                 1647.28340
9               3.38424                 1658.27907
10              3.37221                 1649.01020
11              3.35911                 1639.24568
12              3.35505                 1633.91130
13              3.45891                 1681.02832
14              3.39328                 1645.74032
15              3.39537                 1643.35714
16              3.37342                 1629.35993
17              3.42193                 1649.36833
18              3.43167                 1650.63423
19              3.45611                 1658.93472
20              3.37735                 1617.75161
21              3.35367                 1603.05330
22              3.38395                 1614.14510
23              3.36451                 1601.50724
24              3.36341                 1597.61928
25              3.36488                 1594.95265
26              3.39121                 1604.04470
27              3.44316                 1625.17058
28              3.46199                 1630.59682
29              3.36285                 1580.53809
30              3.38802                 1588.98326
31              3.39398                 1588.38451
32              3.35757                 1567.98332
33              3.35377                 1562.85542
34              3.43302                 1596.35616
35              3.35265                 1555.62867
36              3.37594                 1563.05929
37              3.57434                 1651.34600
38              3.50657                 1616.53015
39              3.47824                 1599.99270
40              3.43989                 1578.90813
41              3.44681                 1578.64127
42              3.36407                 1537.37771
43              3.35832                 1531.39574
44              3.35658                 1527.24208
45              3.35452                 1522.95253
46              3.35567                 1520.11896
47              3.35331                 1515.69838
48              3.35333                 1512.35363
49              3.35215                 1508.46570
50              3.35531                 1506.53329

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_52.png 
51              3.49306                 1564.88864
52              3.35359                 1499.05473
53              3.41925                 1524.98372
54              3.52162                 1567.12090
55              3.49147                 1550.21135
56              3.36783                 1491.94869
57              3.35655                 1483.59687
58              3.40900                 1503.36988
59              3.38128                 1487.76408
60              3.38591                 1486.41581
61              3.35710                 1470.40980
62              3.35437                 1465.85838
63              3.35855                 1464.32954
64              3.35472                 1459.30494
65              3.36141                 1458.85151
66              3.35359                 1452.10360
67              3.38171                 1460.90045
68              3.38606                 1459.39057
69              3.49862                 1504.40660
70              3.38207                 1450.91017
71              3.42677                 1466.65542
72              3.45459                 1475.11078
73              3.35907                 1430.96595
74              3.36526                 1430.23423
75              3.36840                 1428.20330
76              3.36007                 1421.30919
77              3.35401                 1415.39433
78              3.35909                 1414.17563
79              3.35574                 1409.40954
80              3.35900                 1407.42142
81              3.35918                 1404.13849
82              3.35351                 1398.41409
83              3.35346                 1395.03936
84              3.35423                 1392.00587
85              3.36001                 1391.04455
86              3.35387                 1385.15037
87              3.35453                 1382.06595
88              3.37540                 1387.28817
89              3.47493                 1424.72130
90              3.40900                 1394.28264
91              3.35873                 1370.36347
92              3.35239                 1364.42314
93              3.37609                 1370.69295
94              3.44899                 1396.84014
95              3.47202                 1402.69729
96              3.40562                 1372.46607
97              3.38265                 1359.82490
98              3.40264                 1364.45784
99              3.37184                 1348.73720
100             3.37594                 1346.99886

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_54.png 
101             3.57794                 1424.02092
102             3.40611                 1352.22607
103             3.35659                 1329.21083
104             3.37699                 1333.91303
105             3.40049                 1339.79148
106             3.36057                 1320.70558
107             3.39204                 1329.68046
108             3.39272                 1326.55508
109             3.36780                 1313.44200
110             3.37976                 1314.72664
111             3.35748                 1302.70302
112             3.35654                 1298.98214
113             3.40859                 1315.71574
114             3.41997                 1316.68768
115             3.36504                 1292.17344
116             3.35163                 1283.67314
117             3.35371                 1281.11607
118             3.36182                 1280.85494
119             3.35802                 1276.04836
120             3.39025                 1284.90323
121             3.38201                 1278.40129
122             3.38707                 1276.92388
123             3.38986                 1274.58623
124             3.37372                 1265.14575
125             3.42341                 1280.35571
126             3.36864                 1256.50458
127             3.36992                 1253.60987
128             3.44685                 1278.77987
129             3.36274                 1244.21306
130             3.35765                 1238.97285
131             3.36139                 1236.99005
132             3.35780                 1232.31187
133             3.35770                 1228.91893
134             3.35408                 1224.24103
135             3.37260                 1227.62567
136             3.36311                 1220.80784
137             3.36482                 1218.06593
138             3.35968                 1212.84628
139             3.35659                 1208.37312
140             3.44244                 1235.83632
141             3.41490                 1222.53456
142             3.42149                 1221.47193
143             3.45603                 1230.34526
144             3.42956                 1217.49274
145             3.42427                 1212.19229
146             3.45148                 1218.37385
147             3.45342                 1215.60243
148             3.41303                 1197.97493
149             3.43400                 1201.89965
150             3.43859                 1200.06966

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_56.png 
151             3.59457                 1250.91036
152             3.43520                 1192.01544
153             3.41015                 1179.91259
154             3.39632                 1171.73178
155             3.40512                 1171.36128
156             3.43653                 1178.73048
157             3.43215                 1173.79462
158             3.42980                 1169.56282
159             3.43118                 1166.59984
160             3.43850                 1165.65286
161             3.45209                 1166.80676
162             3.38721                 1141.48842
163             3.45100                 1159.53466
164             3.38349                 1133.47015
165             3.42176                 1142.86851
166             3.43391                 1143.49103
167             3.45381                 1146.66359
168             3.44122                 1139.04415
169             3.43990                 1135.16766
170             3.42670                 1127.38266
171             3.41265                 1119.35051
172             3.41608                 1117.05685
173             3.41544                 1113.43474
174             3.41699                 1110.52142
175             3.47908                 1127.22062
176             3.46702                 1119.84875
177             3.44356                 1108.82632
178             3.42794                 1100.36874
179             3.38803                 1084.17024
180             3.43262                 1095.00578
181             3.41756                 1086.78440
182             3.44810                 1093.04833
183             3.43037                 1083.99660
184             3.43997                 1083.59181
185             3.40603                 1069.49342
186             3.37846                 1057.45861
187             3.35036                 1045.31357
188             3.35257                 1042.64927
189             3.35555                 1040.22112
190             3.35072                 1035.37248
191             3.34867                 1031.39005
192             3.35913                 1031.25352
193             3.35308                 1026.04095
194             3.34974                 1021.67223
195             3.34627                 1017.26486
196             3.36614                 1019.94163
197             3.35011                 1011.73292
198             3.34914                 1008.09174
199             3.35072                 1005.21720
200             3.34693                 1000.73058

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_58.png 
Auto-corr tau/N = [0.06144321]
tau/N <= 0.02 = [False]

201             3.49280                 1040.85529
202             3.34585                 993.71864
203             3.34681                 990.65576
204             3.34691                 987.33933
205             3.35233                 985.58502
206             3.35035                 981.65138
207             3.34711                 977.35524
208             3.35587                 976.55730
209             3.34614                 970.38176
210             3.36815                 973.39593
211             3.38941                 976.14950
212             3.35405                 962.61149
213             3.36082                 961.19566
214             3.36996                 960.43945
215             3.37227                 957.72553
216             3.36687                 952.82449
217             3.35640                 946.50508
218             3.35884                 943.83488
219             3.35858                 940.40212
220             3.42666                 956.03898
221             3.50733                 975.03663
222             3.40882                 944.24397
223             3.42843                 946.24668
224             3.41163                 938.19798
225             3.39783                 931.00569
226             3.48435                 951.22673
227             3.45381                 939.43523
228             3.42313                 927.66769
229             3.36870                 909.54765
230             3.39874                 914.26133
231             3.45050                 924.73427
232             3.37358                 900.74639
233             3.47294                 923.80124
234             3.36516                 891.76687
235             3.37203                 890.21513
236             3.41170                 897.27815
237             3.39255                 888.84915
238             3.40962                 889.90978
239             3.41620                 888.21226
240             3.41119                 883.49692
241             3.42002                 882.36464
242             3.37109                 866.36936
243             3.35910                 859.92909
244             3.35958                 856.69341
245             3.36054                 853.57741
246             3.42748                 867.15370
247             3.36069                 846.89438
248             3.36499                 844.61374
249             3.39276                 848.19000
250             3.35910                 836.41565

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_60.png 
Auto-corr tau/N = [0.05382197]
tau/N <= 0.02 = [False]

251             3.52023                 873.01654
252             3.36973                 832.32257
253             3.37223                 829.56981
254             3.36458                 824.32137
255             3.39323                 827.94690
256             3.57015                 867.54621
257             3.58705                 868.06562
258             3.56473                 859.10017
259             3.57104                 857.04888
260             3.45975                 826.88073
261             3.43762                 818.15356
262             3.39460                 804.52091
263             3.38662                 799.24208
264             3.38960                 796.55624
265             3.37064                 788.72906
266             3.45290                 804.52640
267             3.37452                 782.88818
268             3.40876                 787.42264
269             3.47186                 798.52826
270             3.39486                 777.42408
271             3.46769                 790.63378
272             3.51063                 796.91392
273             3.42144                 773.24521
274             3.35905                 755.78512
275             3.35843                 752.28854
276             3.43080                 765.06862
277             3.38239                 750.89102
278             3.39837                 751.03977
279             3.36005                 739.21056
280             3.39065                 742.55213
281             3.53311                 770.21733
282             3.39285                 736.24802
283             3.39414                 733.13359
284             3.47280                 746.65265
285             3.45608                 739.60219
286             3.44084                 732.89828
287             3.37691                 715.90534
288             3.44159                 726.17655
289             3.36503                 706.65567
290             3.40647                 711.95160
291             3.35604                 698.05590
292             3.35623                 694.74023
293             3.35294                 690.70482
294             3.35366                 687.49989
295             3.35369                 684.15215
296             3.35331                 680.72091
297             3.35628                 677.96816
298             3.35766                 674.88906
299             3.35847                 671.69340
300             3.41031                 678.65089

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_62.png 
Auto-corr tau/N = [0.0455079]
tau/N <= 0.02 = [False]

301             3.56518                 705.90643
302             3.38892                 667.61724
303             3.38312                 663.09093
304             3.40532                 664.03720
305             3.44714                 668.74555
306             3.41147                 658.41352
307             3.36019                 645.15648
308             3.35450                 640.70874
309             3.35622                 637.68218
310             3.45529                 653.04981
311             3.43378                 645.55083
312             3.42457                 640.39440
313             3.35281                 623.62266
314             3.35397                 620.48371
315             3.34991                 616.38344
316             3.50427                 641.28104
317             3.61435                 657.81188
318             3.37181                 610.29815
319             3.43519                 618.33456
320             3.45710                 618.82018
321             3.43419                 611.28564
322             3.39766                 601.38547
323             3.35299                 590.12606
324             3.36249                 588.43540
325             3.40951                 593.25439
326             3.35907                 581.11876
327             3.38167                 581.64655
328             3.47247                 593.79203
329             3.40528                 578.89828
330             3.53675                 597.71159
331             3.50958                 589.60944
332             3.48264                 581.60055
333             3.43375                 570.00300
334             3.42031                 564.35049
335             3.43581                 563.47350
336             3.47534                 566.48042
337             3.48137                 563.98145
338             3.54155                 570.18891
339             3.58477                 573.56384
340             3.49435                 555.60213
341             3.46121                 546.87118
342             3.44965                 541.59458
343             3.46415                 540.40724
344             3.40744                 528.15320
345             3.45574                 532.18458
346             3.43384                 525.37829
347             3.47432                 528.09618
348             3.41621                 515.84711
349             3.48263                 522.39405
350             3.47238                 517.38402

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_64.png 
Auto-corr tau/N = [0.04246048]
tau/N <= 0.02 = [False]

351             4.08757                 604.95992
352             3.45484                 507.86119
353             3.43072                 500.88454
354             3.47324                 503.61922
355             3.46237                 498.58157
356             3.49071                 499.17139
357             3.43709                 488.06721
358             3.45032                 486.49498
359             3.41455                 478.03728
360             3.43644                 477.66474
361             3.43086                 473.45937
362             3.47193                 475.65441
363             3.46565                 471.32867
364             3.49823                 472.26038
365             3.47480                 465.62360
366             3.45101                 458.98460
367             3.42686                 452.34486
368             3.42379                 448.51675
369             3.57062                 464.17995
370             3.42973                 442.43581
371             3.52012                 450.57472
372             3.50694                 445.38138
373             3.43249                 432.49424
374             3.50803                 438.50337
375             3.54904                 440.08108
376             3.49957                 430.44736
377             3.45221                 421.16913
378             3.42890                 414.89714
379             3.43260                 411.91164
380             3.37789                 401.96915
381             3.38376                 399.28356
382             3.42072                 400.22412
383             3.49280                 405.16457
384             3.49393                 401.80160
385             3.71224                 423.19513
386             3.64385                 411.75494
387             3.51006                 393.12683
388             3.58802                 398.26978
389             3.49705                 384.67572
390             3.43854                 374.80140
391             3.53080                 381.32618
392             3.46543                 370.80090
393             3.42678                 363.23857
394             3.44966                 362.21462
395             3.58761                 373.11113
396             3.46865                 357.27116
397             3.44685                 351.57839
398             3.50520                 354.02550
399             3.46028                 346.02830
400             3.45796                 342.33844

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_66.png 
Auto-corr tau/N = [0.04080112]
tau/N <= 0.02 = [False]

401             3.57711                 350.55707
402             3.51762                 341.20943
403             3.42937                 329.21990
404             3.46937                 329.59015
405             3.47811                 326.94225
406             3.47332                 323.01867
407             3.44257                 316.71598
408             3.43204                 312.31600
409             3.43595                 309.23550
410             3.45421                 307.42433
411             3.40084                 299.27392
412             3.40612                 296.33288
413             3.42498                 294.54802
414             3.44576                 292.88977
415             3.42702                 287.86976
416             3.39831                 282.05990
417             3.47356                 284.83167
418             3.42902                 277.75070
419             3.41093                 272.87400
420             3.56323                 281.49478
421             3.42079                 266.82139
422             3.47361                 267.46805
423             3.45118                 262.28930
424             3.42843                 257.13263
425             3.46243                 256.22012
426             3.42727                 250.19056
427             3.40093                 244.86710
428             3.40649                 241.86107
429             3.41663                 239.16410
430             3.42946                 236.63308
431             3.42085                 232.61753
432             3.40804                 228.33901
433             3.47319                 229.23054
434             3.43410                 223.21618
435             3.42110                 218.95014
436             3.38678                 213.36714
437             3.40735                 211.25545
438             3.38928                 206.74590
439             3.41785                 205.07076
440             3.47794                 205.19864
441             3.42152                 198.44804
442             3.38243                 192.79874
443             3.43581                 192.40553
444             3.43627                 188.99485
445             3.45305                 186.46454
446             3.47383                 184.11315
447             3.42361                 178.02782
448             3.43519                 175.19489
449             3.43372                 171.68610
450             3.41981                 167.57059

 ac convergence test in progress...
../_images/tutorials_07_ifs_psfsub_fm_planets_215_68.png 
Auto-corr tau/N = [0.03893467]
tau/N <= 0.02 = [False]

451             3.59928                 172.76558
452             3.45187                 162.23794
453             3.42781                 157.67926
454             3.45574                 155.50843
455             3.63760                 160.05422
456             3.40545                 146.43439
457             3.41887                 143.59258
458             3.49040                 143.10628
459             3.45107                 138.04288
460             3.42442                 133.55254
461             3.44353                 130.85410
462             3.38884                 125.38708
463             3.43365                 123.61147
464             3.46346                 121.22114
465             3.42119                 116.32056
466             3.40654                 112.41575
467             3.43290                 109.85293
468             3.42002                 106.02078
469             3.44056                 103.21677
470             3.41010                 98.89284
471             3.43815                 96.26806
472             3.45125                 93.18378
473             3.42340                 89.00850
474             3.41715                 85.42865
475             3.45000                 82.80000
476             3.38471                 77.84835
477             3.41100                 75.04207
478             3.43790                 72.19586
479             3.43830                 68.76592
480             3.41881                 64.95739
481             3.42715                 61.68870
482             3.39495                 57.71413
483             3.41251                 54.60011
484             3.44851                 51.72772
485             3.47022                 48.58312
486             3.40442                 44.25743
487             3.40580                 40.86961
488             3.42750                 37.70247
489             3.44953                 34.49529
490             3.49915                 31.49235
491             3.45939                 27.67508
492             3.45868                 24.21074
493             3.45461                 20.72769
494             3.47143                 17.35717
495             3.43434                 13.73737
496             3.43310                 10.29930
497             3.39262                 6.78525
498             3.37723                 3.37723
499             3.43388                 0.00000
We have reached the limit # of steps without convergence
Running time:  0:28:33.805285
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

[92]:

chain = []
for i, ch in enumerate(test_ch):

    with open('../datasets/MCMC_results_pl{}_ch{}'.format(lab, ch), 'rb') as fi:
        myPickler = pickle.Unpickler(fi)
        mcmc_result = myPickler.load()

    chain.append(mcmc_result['chain'])

Let’s visualize the chain after setting labels for each parameter:

[93]:

for i, ch in enumerate(test_ch):
    labels = [r"$f$ (ch.{:.0f})".format(ch)]
    show_walk_plot(chain[i], labels=labels)
../_images/tutorials_07_ifs_psfsub_fm_planets_218_0.png
../_images/tutorials_07_ifs_psfsub_fm_planets_218_1.png
../_images/tutorials_07_ifs_psfsub_fm_planets_218_2.png

Let’s visualize the corner plot after burning the first 30% of the chain:

[94]:

burnin = 0.3
for i, ch in enumerate(test_ch):
    labels = [r"$f$ (ch.{:.0f})".format(ch)]
    show_corner_plot(chain[i], burnin=burnin, labels=labels)
../_images/tutorials_07_ifs_psfsub_fm_planets_220_0.png
../_images/tutorials_07_ifs_psfsub_fm_planets_220_1.png
../_images/tutorials_07_ifs_psfsub_fm_planets_220_2.png

To calculate 1-sigma (68%) confidence intervals on each parameter, let’s first flatten the chain:

[95]:

npar = 1

[96]:

isamples_flat = chain[0][:, int(chain[0].shape[1]//(1/burnin)):, :].reshape((-1,npar))
for i in range(1,len(test_ch)):
    isamples_flat = np.hstack((isamples_flat,
                               chain[i][:, int(chain[i].shape[1]//(1/burnin)):, :].reshape((-1,npar))))

Then use the confidence function.

Below, the first set of plots show the estimate, ground truth value and 68% confidence interval (in shaded area).

By setting ‘gaussian_fit’ to True, a second set of plots is also shown, where a Gaussian function is fit to each posterior distribution in order to infer mean and standard deviation for each parameter.

Note that the function can accept a ground truth value gt, which will be shown in the plots if provided.

[97]:

gt = []
labels = []
for i, ch in enumerate(test_ch):
    gt.append(flux_c_scal[ch])
    labels.append(r"$f$ (ch.{:.0f})".format(ch))

mu, sigma = confidence(isamples_flat, cfd=68, gaussian_fit=True, verbose=False,
                       gt=gt, save=False, title=True, labels=labels)
../_images/tutorials_07_ifs_psfsub_fm_planets_225_0.png

Overall we see that the estimated fluxes are consistent with the ground truth values (shown with dashed blue lines).