From b4af6e05c6eacfd9f30cd168ac73d6ea515a2cee Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Wed, 11 May 2022 18:35:38 +0200 Subject: [PATCH 01/17] Added automatic tests for NMF-ADI, PCA-RDI (mask), andromeda, and recentering functions --- tests/test_pipeline_adi.py | 4 + tests/test_pipeline_adi_andro.py | 8 ++ tests/test_pipeline_rdi.py | 1 + tests/test_preproc_recentering.py | 144 ++++++++++++++++++++++++------ vip_hci/preproc/recentering.py | 13 ++- 5 files changed, 137 insertions(+), 33 deletions(-) diff --git a/tests/test_pipeline_adi.py b/tests/test_pipeline_adi.py index b3b20763..a98485f6 100644 --- a/tests/test_pipeline_adi.py +++ b/tests/test_pipeline_adi.py @@ -77,6 +77,8 @@ def algo_nmf_annular(ds): def algo_pca(ds): return vip.psfsub.pca(ds.cube, ds.angles, svd_mode='arpack') +def algo_pca_linalg(ds): + return vip.psfsub.pca(ds.cube, ds.angles, svd_mode='eigen') def algo_pca_drot(ds): return vip.psfsub.pca(ds.cube, ds.angles, ncomp=4, fwhm=ds.fwhm, @@ -155,11 +157,13 @@ def verify_expcoord(vectory, vectorx, exp_yx): (algo_medsub_annular, snrmap_fast), (algo_xloci, snrmap_fast), (algo_nmf, snrmap_fast), + (algo_nmf_drot, snrmap_fast), (algo_nmf_annular, snrmap_fast), (algo_llsg, snrmap_fast), (algo_frdiff, snrmap_fast), (algo_frdiff4, snrmap_fast), (algo_pca, snrmap_fast), + (algo_pca_linalg, snrmap_fast), (algo_pca_drot, snrmap_fast), (algo_pca_cevr, snrmap_fast), (algo_pca_grid, snrmap_fast), diff --git a/tests/test_pipeline_adi_andro.py b/tests/test_pipeline_adi_andro.py index 6829352b..9942355c 100644 --- a/tests/test_pipeline_adi_andro.py +++ b/tests/test_pipeline_adi_andro.py @@ -43,6 +43,13 @@ def algo_andromeda(ds): contrast, snr, snr_n, stdcontrast, stdcontrast_n, likelihood, r = res return snr_n +def algo_andromeda_l1(ds): + res = vip.invprob.andromeda(ds.cube[:, :-1, :-1], oversampling_fact=1.8, + angles=ds.angles, psf=ds.psf, iwa=0.5, + opt_method='l1') + contrast, snr, snr_n, stdcontrast, stdcontrast_n, likelihood, r = res + return snr_n + def algo_andromeda_fast(ds): res = vip.invprob.andromeda(ds.cube[:, :-1, :-1], oversampling_fact=0.5, @@ -98,6 +105,7 @@ def verify_expcoord(vectory, vectorx, exp_yx): @parametrize("algo, make_detmap", [ (algo_andromeda, None), + (algo_andromeda_l1, None), (algo_andromeda_fast, None) # (algo_fast_paco, None), # (algo_fast_paco_parallel, None), diff --git a/tests/test_pipeline_rdi.py b/tests/test_pipeline_rdi.py index 160260b7..f2059ea5 100644 --- a/tests/test_pipeline_rdi.py +++ b/tests/test_pipeline_rdi.py @@ -121,6 +121,7 @@ def verify_expcoord(vectory, vectorx, exp_yx): @parametrize("algo, make_detmap", [ (algo_pca, snrmap_fast), + (algo_pca_mask, snrmap_fast), (algo_pca_annular_ardi, snrmap_fast), (algo_pca_annular_rdi, snrmap_fast), (algo_nmf, snrmap_fast), diff --git a/tests/test_preproc_recentering.py b/tests/test_preproc_recentering.py index b32d81ef..884caf4b 100644 --- a/tests/test_preproc_recentering.py +++ b/tests/test_preproc_recentering.py @@ -1,10 +1,13 @@ from sklearn.metrics import mean_squared_error from vip_hci.preproc import (cube_recenter_2dfit, cube_recenter_dft_upsampling, cube_recenter_satspots, cube_recenter_via_speckles, - frame_shift) + cube_recenter_radon, frame_shift, cube_subsample, + cube_correct_nan, approx_stellar_position) +from vip_hci.var import frame_center import vip_hci as vip import hciplot from astropy.modeling import models +from astropy.utils.data import download_file from .helpers import fixture, np import copy import matplotlib as mpl @@ -66,6 +69,25 @@ def get_cube(example_dataset_adi): return dsi +@fixture(scope="module") +def get_ifs_cube(example_dataset_sdi): + """ + Get the ADI+IFS sequence from conftest.py. + + Parameters + ---------- + example_dataset_sdi : fixture + Taken automatically from ``conftest.py``. + + Returns + ------- + dsi : VIP Dataset + + """ + dsi = copy.copy(example_dataset_sdi) + + return dsi + def shift_cube(cube, randax, randay): return np.array([frame_shift(cube[i], randay[i], randax[i]) @@ -198,8 +220,35 @@ def do_recenter(method, cube, shiftx, shifty, errormsg, mse=1e-2, # d88P"Y88b Y88b. Y8b. X88 Y88b. X88 # dP" "Yb "Y888 "Y8888 88888P' "Y888 88888P' +def test_approx_star(debug=True): + """ + tests `approx_stellar_position`. The data cube is generated from a 2D + gaussian (positive case). + """ + global seed + + if debug: + html("

===== test_2d =====

") + + method = approx_stellar_position + errormsg = 'Error when recentering with approx. stellar position method' + n_frames = 6 + + + for model, name in {"gauss": "Gaussian"}.items(): -def test_2d(debug=False): + # ===== odd + cube = create_cube_with_gauss2d(shape=(n_frames, 9, 9), mean=4, + stddev=1) + cube += np.random.normal(0, 0.1, cube.shape) + cy, cx = frame_center(cube) + esty, estx = method(cube, fwhm=1) + mse = 0.02 + assert mean_squared_error([cy]*len(esty), esty) < mse, errormsg + assert mean_squared_error([cx]*len(estx), estx) < mse, errormsg + + +def test_2d(debug=True): """ tests `cube_recenter_2dfit`. The data cube is generated from a 2D gaussian (positive case) or a 2D "ring" (difference of two 2D gaussians, negative @@ -223,8 +272,8 @@ def test_2d(debug=False): cube = create_cube_with_gauss2d(shape=(n_frames, 9, 9), mean=4, stddev=1) - method_args = dict(fwhm=1, subi_size=5, model=model, verbose=False, - negative=False, full_output=True, plot=False) + method_args = dict(fwhm=1, subi_size=5, model=model, verbose=True, + negative=False, full_output=True, plot=True) do_recenter(method, cube, randax, randay, errormsg=errormsg.format(name), debug=debug, **method_args) @@ -232,8 +281,8 @@ def test_2d(debug=False): cube = create_cube_with_gauss2d(shape=(n_frames, 10, 10), mean=5, stddev=1) - method_args = dict(fwhm=1, subi_size=6, model=model, verbose=False, - negative=False, full_output=True, plot=False) + method_args = dict(fwhm=1, subi_size=6, model=model, verbose=True, + negative=False, full_output=True, plot=True) do_recenter(method, cube, randax, randay, errormsg=errormsg.format(name), debug=debug, **method_args) @@ -241,8 +290,8 @@ def test_2d(debug=False): cube = create_cube_with_gauss2d_ring(shape=(n_frames, 9, 9), mean=4, stddev_outer=3, stddev_inner=2) - method_args = dict(fwhm=1, subi_size=5, model=model, verbose=False, - negative=True, full_output=True, plot=False) + method_args = dict(fwhm=1, subi_size=5, model=model, verbose=True, + negative=True, full_output=True, plot=True) do_recenter(method, cube, randax, randay, errormsg=errormsg.format(name), debug=debug, **method_args) @@ -250,19 +299,19 @@ def test_2d(debug=False): cube = create_cube_with_gauss2d_ring(shape=(n_frames, 10, 10), mean=5, stddev_outer=3, stddev_inner=2) - method_args = dict(fwhm=1, subi_size=6, model=model, verbose=False, - negative=True, full_output=True, plot=False) + method_args = dict(fwhm=1, subi_size=6, model=model, verbose=True, + negative=True, full_output=True, plot=True) do_recenter(method, cube, randax, randay, errormsg=errormsg.format(name), debug=debug, **method_args) - - -def test_dft(debug=False): + + +def test_dft(debug=True): global seed if debug: html("

===== test_dft =====

") method = cube_recenter_dft_upsampling - method_args_additional = dict(verbose=True, full_output=True, plot=False) + method_args_additional = dict(verbose=True, full_output=True, plot=True) errormsg = 'Error when recentering with DFT upsampling method' n_frames = 6 @@ -319,7 +368,7 @@ def test_dft(debug=False): mse_skip_first=True, debug=debug, **method_args) -def test_dft_image(debug=False): +def test_dft_image(debug=True): """ notes: ====== @@ -346,12 +395,12 @@ def test_dft_image(debug=False): # ===== recenter method_args = dict(center_fr1=(51, 51), subi_size=None, verbose=True, - negative=True, full_output=True, plot=False) + negative=True, full_output=True, plot=True) do_recenter(method, cube, randax, randay, errormsg=errormsg, mse_skip_first=True, debug=debug, **method_args) -def test_satspots_image(debug=False): +def test_satspots_image(debug=True): global seed if debug: html("

===== test_satspots_image =====

") @@ -370,13 +419,13 @@ def test_satspots_image(debug=False): # ===== recenter spotcoords = [(41, 109), (109, 109), (41, 41), (109, 41)] # NW NE SW SE - method_args = dict(xy=spotcoords, subi_size=25, plot=False, - full_output=True, verbose=False) + method_args = dict(xy=spotcoords, subi_size=25, mask_center=40, + full_output=True, verbose=True) do_recenter(method, cube, randax, randay, errormsg=errormsg, debug=debug, **method_args) -def test_satspots(debug=False): +def test_satspots(debug=True): global seed if debug: html("

===== test_satspots =====

") @@ -394,16 +443,50 @@ def test_satspots(debug=False): randay = seed.uniform(0, shift_magnitude, size=n_frames) # ===== recenter - method_args = dict(xy=spotcoords, subi_size=9, plot=False, - full_output=True, verbose=False) + method_args = dict(xy=spotcoords, subi_size=9, plot=True, + full_output=True, verbose=True) do_recenter(method, cube, randax, randay, errormsg=errormsg, debug=debug, **method_args) -def test_speckle_recentering(get_cube, debug=False): +def test_radon(debug=True): global seed if debug: - html("

===== test_satspots =====

") + html("

===== test_radon =====

") + + method = cube_recenter_radon + errormsg = 'Error when recentering with Radon transform' + n_frames = 2 + + # ===== datacube + url_prefix = "https://github.com/vortex-exoplanet/VIP_extras/raw/master/datasets" + + f1 = download_file("{}/sphere_ifs_PDS70_cen.fits".format(url_prefix), + cache=True) + + # load fits + cube = vip.fits.open_fits(f1) + + # subsample and correct for NaNs + cube = cube_subsample(cube, 19) + cube = cube_correct_nan(cube) + + # ===== shift + shift_magnitude = 1 + randax = seed.uniform(-shift_magnitude, 0, size=n_frames) + randay = seed.uniform(0, shift_magnitude, size=n_frames) + + # ===== recenter + method_args = dict(hsize=1.0, step=0.1, plot=True, full_output=True, + verbose=True) + do_recenter(method, cube, randax, randay, errormsg=errormsg, debug=debug, + mse=0.04, **method_args) + + +def test_speckle_recentering(get_cube, debug=True): + global seed + if debug: + html("

===== test_speckle_recentering =====

") method = cube_recenter_via_speckles errormsg = 'Error when recentering via speckles' @@ -417,8 +500,11 @@ def test_speckle_recentering(get_cube, debug=False): randay = np.ones(n_frames) # ===== recenter - method_args = dict(plot=False, full_output=True, fwhm=4.2, - recenter_median=True, subframesize=49, imlib='opencv', - interpolation='lanczos4') - do_recenter(method, ds.cube, randax, randay, errormsg=errormsg, debug=debug, - mse=0.2*n_frames, **method_args) + types = ['gaus', 'ann'] + + for ty in types: + method_args = dict(plot=False, full_output=True, fwhm=4.2, fit_type=ty, + recenter_median=True, subframesize=49, imlib='opencv', + interpolation='lanczos4') + do_recenter(method, ds.cube, randax, randay, errormsg=errormsg, + debug=debug, mse=0.04, **method_args) diff --git a/vip_hci/preproc/recentering.py b/vip_hci/preproc/recentering.py index d0189af7..7fc69cbf 100644 --- a/vip_hci/preproc/recentering.py +++ b/vip_hci/preproc/recentering.py @@ -867,7 +867,8 @@ def cube_recenter_radon(array, full_output=False, verbose=True, imlib='vip-fft', cy, cx = frame_center(array[0]) array_rec = array.copy() - for i in Progressbar(range(n_frames), desc="frames", verbose=verbose): + for i in Progressbar(range(n_frames), desc="Recentering frames...", + verbose=verbose): y[i], x[i] = frame_center_radon(array[i], verbose=False, plot=False, **kwargs) array_rec[i] = frame_shift(array[i], cy-y[i], cx-x[i], imlib=imlib, @@ -878,7 +879,7 @@ def cube_recenter_radon(array, full_output=False, verbose=True, imlib='vip-fft', timing(start_time) if full_output: - return array_rec, y, x + return array_rec, cy-y, cx-x else: return array_rec @@ -1541,7 +1542,7 @@ def cube_recenter_via_speckles(cube_sci, cube_ref=None, alignment_iter=5, if recenter_median and fit_type not in {'gaus', 'ann'}: raise TypeError("fit type not recognized. Should be 'ann' or 'gaus'") - if crop and not subframesize < y/2.: + if crop and not subframesize < y: raise ValueError('`Subframesize` is too large') if cube_ref is not None: @@ -1618,7 +1619,11 @@ def cube_recenter_via_speckles(cube_sci, cube_ref=None, alignment_iter=5, else: crop_sz = int(6*fwhm) if not crop_sz % 2: - crop_sz += 1 + # size should be odd and small, but at least 7 for 2D fit + if crop_sz>7: + crop_sz -= 1 + else: + crop_sz += 1 sub_image, y1, x1 = get_square(alignment_cube[0], size=crop_sz, y=ceny, x=cenx, position=True) From b40eae4b5786a6bd15398ab1f97df491d893eabf Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Wed, 11 May 2022 19:27:51 +0200 Subject: [PATCH 02/17] Updated default codecov parameters in yml file --- .codecov.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.codecov.yml b/.codecov.yml index 419985e6..03240796 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -1,6 +1,9 @@ coverage: status: patch: false + range: 50..70 + round: nearest + precision: 2 ignore: - tests/* From d2a402e6c038495780aacf3e7fa02495529fb53f Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Wed, 11 May 2022 22:46:34 +0200 Subject: [PATCH 03/17] tests: minor bug fix in full_output full-frame NMF and turned off debug option in recentering tests --- tests/test_preproc_recentering.py | 16 ++++++++-------- vip_hci/psfsub/nmf_fullfr.py | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/test_preproc_recentering.py b/tests/test_preproc_recentering.py index 884caf4b..0cb93a3c 100644 --- a/tests/test_preproc_recentering.py +++ b/tests/test_preproc_recentering.py @@ -220,7 +220,7 @@ def do_recenter(method, cube, shiftx, shifty, errormsg, mse=1e-2, # d88P"Y88b Y88b. Y8b. X88 Y88b. X88 # dP" "Yb "Y888 "Y8888 88888P' "Y888 88888P' -def test_approx_star(debug=True): +def test_approx_star(debug=False): """ tests `approx_stellar_position`. The data cube is generated from a 2D gaussian (positive case). @@ -248,7 +248,7 @@ def test_approx_star(debug=True): assert mean_squared_error([cx]*len(estx), estx) < mse, errormsg -def test_2d(debug=True): +def test_2d(debug=False): """ tests `cube_recenter_2dfit`. The data cube is generated from a 2D gaussian (positive case) or a 2D "ring" (difference of two 2D gaussians, negative @@ -305,7 +305,7 @@ def test_2d(debug=True): errormsg=errormsg.format(name), debug=debug, **method_args) -def test_dft(debug=True): +def test_dft(debug=False): global seed if debug: html("

===== test_dft =====

") @@ -368,7 +368,7 @@ def test_dft(debug=True): mse_skip_first=True, debug=debug, **method_args) -def test_dft_image(debug=True): +def test_dft_image(debug=False): """ notes: ====== @@ -400,7 +400,7 @@ def test_dft_image(debug=True): mse_skip_first=True, debug=debug, **method_args) -def test_satspots_image(debug=True): +def test_satspots_image(debug=False): global seed if debug: html("

===== test_satspots_image =====

") @@ -425,7 +425,7 @@ def test_satspots_image(debug=True): **method_args) -def test_satspots(debug=True): +def test_satspots(debug=False): global seed if debug: html("

===== test_satspots =====

") @@ -449,7 +449,7 @@ def test_satspots(debug=True): **method_args) -def test_radon(debug=True): +def test_radon(debug=False): global seed if debug: html("

===== test_radon =====

") @@ -483,7 +483,7 @@ def test_radon(debug=True): mse=0.04, **method_args) -def test_speckle_recentering(get_cube, debug=True): +def test_speckle_recentering(get_cube, debug=False): global seed if debug: html("

===== test_speckle_recentering =====

") diff --git a/vip_hci/psfsub/nmf_fullfr.py b/vip_hci/psfsub/nmf_fullfr.py index 0d3d4697..f240a483 100644 --- a/vip_hci/psfsub/nmf_fullfr.py +++ b/vip_hci/psfsub/nmf_fullfr.py @@ -227,13 +227,13 @@ def nmf(cube, angle_list, cube_ref=None, ncomp=1, scaling=None, max_iter=10000, residuals = res_result if handle_neg == 'mask': residuals_cube[fr][yy, xx] = residuals - if fr == n-1: + if fr == n-1 and full_output: for pp in range(ncomp): H_tmp[pp][yy, xx] = H[pp] H = H_tmp else: residuals_cube[fr] = residuals.reshape((y, x)) - if fr == n-1: + if fr == n-1 and full_output: H = H.reshape(ncomp, y, x) if verbose: From 44920b4cae0386d2f920624b4d546ebf2d6f4006 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 12 May 2022 10:22:54 +0200 Subject: [PATCH 04/17] Tests: Minor bug fixes for additional tests on recentering functions --- tests/test_preproc_recentering.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/test_preproc_recentering.py b/tests/test_preproc_recentering.py index 0cb93a3c..301b2f30 100644 --- a/tests/test_preproc_recentering.py +++ b/tests/test_preproc_recentering.py @@ -242,7 +242,9 @@ def test_approx_star(debug=False): stddev=1) cube += np.random.normal(0, 0.1, cube.shape) cy, cx = frame_center(cube) - esty, estx = method(cube, fwhm=1) + est_yx = method(cube, fwhm=1) + esty = est_yx[:,0] + estx = est_yx[:,1] mse = 0.02 assert mean_squared_error([cy]*len(esty), esty) < mse, errormsg assert mean_squared_error([cx]*len(estx), estx) < mse, errormsg @@ -419,7 +421,7 @@ def test_satspots_image(debug=False): # ===== recenter spotcoords = [(41, 109), (109, 109), (41, 41), (109, 41)] # NW NE SW SE - method_args = dict(xy=spotcoords, subi_size=25, mask_center=40, + method_args = dict(xy=spotcoords, subi_size=25, full_output=True, verbose=True) do_recenter(method, cube, randax, randay, errormsg=errormsg, debug=debug, **method_args) @@ -477,7 +479,7 @@ def test_radon(debug=False): randay = seed.uniform(0, shift_magnitude, size=n_frames) # ===== recenter - method_args = dict(hsize=1.0, step=0.1, plot=True, full_output=True, + method_args = dict(hsize=1.0, step=0.1, full_output=True, mask_center=40, verbose=True) do_recenter(method, cube, randax, randay, errormsg=errormsg, debug=debug, mse=0.04, **method_args) From 3c254d03a83939f1bb8d23d916f9e4e59afb79be Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 12 May 2022 12:23:39 +0200 Subject: [PATCH 05/17] Added verbose outputs for completeness curve algorithm --- vip_hci/metrics/completeness.py | 100 +++++++++++++++++++++----------- 1 file changed, 65 insertions(+), 35 deletions(-) diff --git a/vip_hci/metrics/completeness.py b/vip_hci/metrics/completeness.py index a81651ce..c57e5e4f 100644 --- a/vip_hci/metrics/completeness.py +++ b/vip_hci/metrics/completeness.py @@ -4,7 +4,7 @@ Module with completeness curve and map generation function. """ -__author__ = 'C.H. Dahlqvist' +__author__ = 'C.H. Dahlqvist, V. Christiaens' __all__ = ['completeness_curve', 'completeness_map'] @@ -124,13 +124,17 @@ def _estimate_snr_fc(a, b, level, n_fc, cube, psf, angle_list, fwhm, algo, snrmap_fin[indc[0], indc[1]] = 0 max_map = np.nan_to_num(snrmap_fin).max() + if b==2 and max_target-max_map<0: + from hciplot import plot_frames + plot_frames((snrmap_empty, snrmap_temp, snrmap_fin)) + return max_target-max_map, b def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, ini_contrast=None, starphot=1, pxscale=0.1, n_fc=20, completeness=0.95, snr_approximation=True, max_iter=50, - nproc=1, algo_dict={'ncomp': 20}, plot=True, + nproc=1, algo_dict={}, verbose=True, plot=True, dpi=vip_figdpi, save_plot=None, object_name=None, fix_y_lim=(), figsize=vip_figsize): """ @@ -204,6 +208,8 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, algo_dict: dictionary, optional Any other valid parameter of the post-processing algorithms can be passed here, including e.g. imlib and interpolation. + verbose : bool, optional + Whether to print more info while running the algorithm. Default: True. plot : bool, optional Whether to plot the final contrast curve or not. True by default. dpi : int optional @@ -275,6 +281,8 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, idx = find_nearest(ini_rads, ad) ini_contrast.append(ini_cc[idx]) + if verbose: + print("Calculating initial SNR map with no injected companion...") argl = getfullargspec(algo).args if 'cube' in argl and 'angle_list' in argl: if 'fwhm' in argl: @@ -303,19 +311,22 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, psf.shape[1])) for k in range(len(an_dist)): - a = an_dist[k] level = ini_contrast[k] pos_detect = [] + if verbose: + print("*** Calculating contrast at r = {} ***".format(a)) + detect_bound = [None, None] level_bound = [None, None] - it = 0 + ii = 0 err_msg = "Could not converge on a contrast level matching required " err_msg += "completeness within {} iterations. Tested level: {}. " - err_msg += "Consider increasing min radius." + err_msg += "Is there too much self-subtraction? Consider decreasing " + err_msg += "ncomp if using PCA, or increasing minimum requested radius." - while len(pos_detect) == 0 and it < max_iter: + while len(pos_detect) == 0 and ii < max_iter: pos_detect = [] pos_non_detect = [] val_detect = [] @@ -323,7 +334,8 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, res = pool_map(nproc, _estimate_snr_fc, a, iterable(range(0, n_fc)), level, n_fc, cube, psf, angle_list, fwhm, algo, - algo_dict, snrmap_empty, starphot, approximated=True) + algo_dict, snrmap_empty, starphot, + approximated=snr_approximation) for res_i in res: @@ -336,9 +348,13 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, if len(pos_detect) == 0: level = level*1.5 - it += 1 + ii += 1 - if it == max_iter: + if verbose: + msg = "Found contrast level for first TP detection: {}" + print(msg.format(level)) + + if ii == max_iter: raise ValueError(err_msg.format(max_iter, level)) if len(pos_detect) > round(completeness*n_fc): @@ -355,8 +371,8 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, cond1 = (detect_bound[0] is None or detect_bound[1] is None) cond2 = (len(pos_detect) != round(completeness*n_fc)) - it = 0 - while cond1 and cond2 and it 0: pos_detect.append(res_i[1]) val_detect.append(res_i[0]) @@ -425,9 +440,13 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, cond1 = (detect_bound[0] is None or detect_bound[1] is None) cond2 = (len(pos_detect) != comp_temp) - it += 1 + ii += 1 + + if verbose: + msg = "Found lower and upper bounds of sought contrast: {}" + print(msg.format(level_bound)) - if it == max_iter: + if ii == max_iter: raise ValueError(err_msg.format(max_iter, level)) if len(pos_detect) != round(completeness*n_fc): @@ -437,8 +456,8 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, pos_detect = pos_detect_temp.copy() val_detect = val_detect_temp.copy() - it = 0 - while len(pos_detect) != round(completeness*n_fc) and it < max_iter: + ii = 0 + while len(pos_detect) != round(completeness*n_fc) and ii < max_iter: fact = (level_bound[1]-level_bound[0]) / \ (detect_bound[1]-detect_bound[0]) level = level_bound[0]+fact*(completeness*n_fc-detect_bound[0]) @@ -446,7 +465,8 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, res = pool_map(nproc, _estimate_snr_fc, a, iterable(-np.sort(-np.array(pos_non_detect))), level, n_fc, cube, psf, angle_list, fwhm, algo, - algo_dict, snrmap_empty, starphot, approximated=True) + algo_dict, snrmap_empty, starphot, + approximated=snr_approximation) it = len(pos_non_detect)-1 for res_i in res: @@ -476,13 +496,14 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, val_non_detect = val_non_detect_temp.copy() pos_detect = pos_detect_temp.copy() val_detect = val_detect_temp.copy() - it += 1 + ii += 1 - if it == max_iter: + if ii == max_iter: raise ValueError(err_msg.format(max_iter, level)) - - print("Distance: "+"{}".format(a)+" Final contrast " + - "{}".format(level)) + + if verbose: + msg = "=> found final contrast for {}% completeness: {}" + print(msg.format(completeness*100, level)) cont_curve[k] = level # plotting @@ -533,7 +554,7 @@ def completeness_curve(cube, angle_list, psf, fwhm, algo, an_dist=None, def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, starphot=1, n_fc=20, snr_approximation=True, nproc=1, - algo_dict={'ncomp': 20}, verbose=True): + algo_dict={}, verbose=True): """ Function allowing the computation of two dimensional (radius and completeness) contrast curves with any psf-subtraction algorithm provided by @@ -592,7 +613,7 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, as 1/n_fc and the maximum is 1-1/n_fc). Default 20. snr_approximated : bool, optional If True, an approximated S/N map is generated. If False the - approach of Mawett et al. is used (2014). Default is True + approach of Mawet et al. is used (2014). Default is True nproc : int or None Number of processes for parallel computing. algo_dict @@ -606,9 +627,12 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, ------- an_dist: 1d numpy ndarray Radial distances where the contrasts are calculated + comp_levels: 1d numpy ndarray + Completeness levels for which the contrasts are calculated cont_curve: 2d numpy ndarray Contrast matrix, with the first axis associated to the radial distances - and the second axis associated to the completeness level. + and the second axis associated to the completeness level, calculated + from 1/n_fc to (n_fc-1)/n_fc. """ @@ -688,7 +712,8 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, pos_non_detect = [] res = pool_map(nproc, _estimate_snr_fc, a, iterable(range(0, n_fc)), level, n_fc, cube, psf, angle_list, fwhm, algo, - algo_dict, snrmap_empty, starphot, approximated=True) + algo_dict, snrmap_empty, starphot, + approximated=snr_approximation) for res_i in res: @@ -709,7 +734,8 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, res = pool_map(nproc, _estimate_snr_fc, a, iterable(-np.sort(-np.array(pos_detect))), level, n_fc, cube, psf, angle_list, fwhm, algo, algo_dict, - snrmap_empty, starphot, approximated=True) + snrmap_empty, starphot, + approximated=snr_approximation) it = len(pos_detect)-1 for res_i in res: @@ -724,7 +750,7 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, list(pos_non_detect.copy())] if verbose: - print("Lower boundary found") + print("Lower boundary ({:.0f}%) found: {}".format(100/n_fc, level)) level = contrast_matrix[k, np.where(contrast_matrix[k, :] > 0)[0][-1]] @@ -738,7 +764,8 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, res = pool_map(nproc, _estimate_snr_fc, a, iterable(-np.sort(-np.array(pos_non_detect))), level, n_fc, cube, psf, angle_list, fwhm, algo, - algo_dict, snrmap_empty, starphot, approximated=True) + algo_dict, snrmap_empty, starphot, + approximated=snr_approximation) it = len(pos_non_detect)-1 for res_i in res: @@ -753,7 +780,8 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, list(pos_non_detect.copy())] if verbose: - print("Upper boundary found") + print("Upper boundary ({:.0f}%) found: {}".format(100*(n_fc-1)/n_fc, + level)) missing = np.where(contrast_matrix[k, :] == 0)[0] computed = np.where(contrast_matrix[k, :] > 0)[0] @@ -784,7 +812,7 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, iterable(-np.sort(-np.array(pos_detect))), level, n_fc, cube, psf, angle_list, fwhm, algo, algo_dict, snrmap_empty, starphot, - approximated=True) + approximated=snr_approximation) it = len(pos_detect)-1 for res_i in res: @@ -808,7 +836,7 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, iterable(-np.sort(-np.array(pos_non_detect))), level, n_fc, cube, psf, angle_list, fwhm, algo, algo_dict, snrmap_empty, starphot, - approximated=True) + approximated=snr_approximation) it = len(pos_non_detect)-1 for res_i in res: @@ -838,4 +866,6 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, computed = np.where(contrast_matrix[k, :] > 0)[0] missing = np.where(contrast_matrix[k, :] == 0)[0] - return an_dist, contrast_matrix + comp_levels = np.linspace(1/n_fc, 1-1/n_fc, n_fc-1, endpoint=True) + + return an_dist, comp_levels, contrast_matrix[1:-1] From db25fbeb3058d46eb93a2ee4464a38d68d67162c Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 12 May 2022 12:25:43 +0200 Subject: [PATCH 06/17] typo fix in completeness curve --- vip_hci/metrics/completeness.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vip_hci/metrics/completeness.py b/vip_hci/metrics/completeness.py index c57e5e4f..61f9be46 100644 --- a/vip_hci/metrics/completeness.py +++ b/vip_hci/metrics/completeness.py @@ -750,7 +750,7 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, list(pos_non_detect.copy())] if verbose: - print("Lower boundary ({:.0f}%) found: {}".format(100/n_fc, level)) + print("Lower bound ({:.0f}%) found: {}".format(100/n_fc, level)) level = contrast_matrix[k, np.where(contrast_matrix[k, :] > 0)[0][-1]] @@ -780,8 +780,8 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, list(pos_non_detect.copy())] if verbose: - print("Upper boundary ({:.0f}%) found: {}".format(100*(n_fc-1)/n_fc, - level)) + print("Upper bound ({:.0f}%) found: {}".format(100*(n_fc-1)/n_fc, + level)) missing = np.where(contrast_matrix[k, :] == 0)[0] computed = np.where(contrast_matrix[k, :] > 0)[0] From acfe9ca8bba7929e45d1a77ab9589e565cc7ffa2 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 12 May 2022 12:28:58 +0200 Subject: [PATCH 07/17] fixed contrast matrix output format to exclude 0 and 1 completeness --- vip_hci/metrics/completeness.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vip_hci/metrics/completeness.py b/vip_hci/metrics/completeness.py index 61f9be46..0064a9fb 100644 --- a/vip_hci/metrics/completeness.py +++ b/vip_hci/metrics/completeness.py @@ -868,4 +868,4 @@ def completeness_map(cube, angle_list, psf, fwhm, algo, an_dist, ini_contrast, comp_levels = np.linspace(1/n_fc, 1-1/n_fc, n_fc-1, endpoint=True) - return an_dist, comp_levels, contrast_matrix[1:-1] + return an_dist, comp_levels, contrast_matrix[:, 1:-1] From 7458bdb8fbcc870a24da0c3d5705e09216415f65 Mon Sep 17 00:00:00 2001 From: VChristiaens Date: Thu, 12 May 2022 12:37:46 +0200 Subject: [PATCH 08/17] tuto4: updated completeness curves and maps example --- docs/source/tutorials/04_metrics.ipynb | 183 ++++++++++++++----------- 1 file changed, 102 insertions(+), 81 deletions(-) diff --git a/docs/source/tutorials/04_metrics.ipynb b/docs/source/tutorials/04_metrics.ipynb index 38878ac6..d3d7de15 100644 --- a/docs/source/tutorials/04_metrics.ipynb +++ b/docs/source/tutorials/04_metrics.ipynb @@ -13,7 +13,7 @@ "source": [ "> Authors: *Carlos Alberto Gomez Gonzalez* and *Valentin Christiaens* \n", "> Suitable for VIP *v1.0.0* onwards \n", - "> Last update: *2022/05/04*" + "> Last update: *2022/05/12*" ] }, { @@ -2376,12 +2376,12 @@ "\n", "\n", "\n", - "
\n", + "
\n", "\n", "