From 74ccfe8cbc5406d6c64586e37afb9e6f2f0d7e7f Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Tue, 13 Dec 2022 11:04:35 +0800 Subject: [PATCH 01/12] cube_correct_nan updated with multiprocessing --- vip_hci/preproc/cosmetics.py | 43 +++++++++++++++++++++++++++++++----- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/vip_hci/preproc/cosmetics.py b/vip_hci/preproc/cosmetics.py index 07fad97d..f454bf12 100644 --- a/vip_hci/preproc/cosmetics.py +++ b/vip_hci/preproc/cosmetics.py @@ -1,5 +1,5 @@ #! /usr/bin/env python - +# check123 """ Module with cosmetics procedures. Contains the function for bad pixel fixing. Also functions for cropping cubes. @@ -21,6 +21,9 @@ from ..config.utils_conf import pool_map, iterable from ..stats import sigma_filter from ..var import frame_center, get_square +from multiprocessing import Process +import multiprocessing +from multiprocessing import shared_memory, set_start_method def cube_crop_frames(array, size, xy=None, force=False, verbose=True, @@ -346,11 +349,41 @@ def cube_correct_nan(cube, neighbor_box=3, min_neighbors=3, verbose=False, print(msg.format(zz, nnanpix)) else: if verbose: - msg = "Correcting NaNs in multiprocessing..." + msg = "Correcting NaNs in multiprocessing using ADACS' approach..." print(msg) - res = pool_map(nproc, nan_corr_2d, iterable(obj_tmp), neighbor_box, - min_neighbors, half_res_y, verbose, False) - obj_tmp = np.array(res, dtype=object) + #creation of shared memory blob + shm = shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) + #creating an array object similar to obj_tmp and occupying shared memory. + obj_tmp_shared=np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm.buf) + #nan_corr_2d_mp function passes the frame and other details to nan_corr_2d for processing. + def nan_corr_2d_mp(i,obj, neighbor_box, min_neighbors, half_res_y, verbose): + obj_tmp_shared[i]=nan_corr_2d(obj, neighbor_box, min_neighbors, half_res_y,verbose,full_output=False) + # _nan_corr_2d unfolds the arguments + global _nan_corr_2d_mp + def _nan_corr_2d_mp(args): + nan_corr_2d_mp(*args) + + context=multiprocessing.get_context('fork') + #processes argumnet is not provided to pool as no. of processes is + #inferred from os.cpu_count() + pool=context.Pool(processes=nproc, maxtasksperchild=1) + #creating a list of arguments + args=[] + for j in range(n_z): + args.append([j, obj_tmp[j], neighbor_box, min_neighbors, half_res_y, verbose]) + try: + pool.map_async(_nan_corr_2d_mp,args, chunksize=1 ).get(timeout=10_000_000) + finally: + pool.close() + pool.join() + obj_tmp[:]=obj_tmp_shared[:] + shm.close() + shm.unlink() + #Original Multiprocessing code commented out below and ADACS approach is activated. + # res = pool_map(nproc, nan_corr_2d, iterable(obj_tmp), neighbor_box, + # min_neighbors, half_res_y, verbose, False) + # obj_tmp = np.array(res, dtype=object) + if verbose: print('All nan pixels are corrected.') From c319f33788d9dcc3b26d0bef50e2624d1fc84e6d Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Tue, 13 Dec 2022 17:44:26 +0800 Subject: [PATCH 02/12] multiprocessing added on cube_fix_badpix_clump --- vip_hci/preproc/badpixremoval.py | 120 +++++++++++++++++++++++++------ 1 file changed, 100 insertions(+), 20 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 874d3d0f..bd70bfc8 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -34,6 +34,9 @@ from ..config.utils_conf import pool_map, iterable from .rescaling import find_scal_vector, frame_rescaling from .cosmetics import frame_pad +from multiprocessing import Process +import multiprocessing +from multiprocessing import shared_memory, set_start_method import warnings try: @@ -592,7 +595,7 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, r_in_std, def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, cx=None, fwhm=4., sig=4., protect_mask=0, half_res_y=False, min_thr=None, max_nit=15, mad=True, - verbose=True, full_output=False): + verbose=True, full_output=False, nproc=1): """ Function to identify and correct clumps of bad pixels. Very fast when a bad pixel map is provided. If a bad pixel map is not provided, the bad pixel @@ -656,6 +659,10 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, full_output: bool, {False,True}, optional Whether to return as well the cube of bad pixel maps and the cube of defined annuli. + nproc: int, optional + This feature is added following ADACS update. Refers to the number of processors + available for calculations. Choosing a number >1 enables multiprocessing for the + correction of frames. Returns ------- @@ -806,16 +813,54 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, cx = [cx]*n_z if isinstance(fwhm, (float, int)): fwhm = [fwhm]*n_z - bpix_map_cumul = np.zeros_like(obj_tmp) - for i in range(n_z): - if verbose: - print('************Frame # ', i, ' *************') - obj_tmp[i], bpix_map_cumul[i] = bp_removal_2d(obj_tmp[i], cy[i], - cx[i], fwhm[i], - sig, protect_mask, - bpm_mask, min_thr, - half_res_y, mad, - verbose) + if nproc==1: + bpix_map_cumul = np.zeros_like(obj_tmp) + for i in range(n_z): + if verbose: + print('************Frame # ', i, ' *************') + obj_tmp[i], bpix_map_cumul[i] = bp_removal_2d(obj_tmp[i], cy[i], + cx[i], fwhm[i], + sig, protect_mask, + bpm_mask, min_thr, + half_res_y, mad, + verbose) + else: + msg="Cleaning frames using ADACS' multiprocessing appraoch" + print(msg) + #creating shared memory buffer space for the image cube. + shm_clump= shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) + obj_tmp_shared_clump = np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump.buf) + #creating shared memory buffer space for the bad pixel cube. + shm_clump_bpix= shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) + #works with dtype=obj_tmp.dtype but not dtype=int + bpix_map_cumul_shared= np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump_bpix.buf) + + def mp_clump_slow(j, obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose): + obj_tmp_shared_clump[j], bpix_map_cumul_shared[j] = bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose) + + global _mp_clump_slow + + def _mp_clump_slow(args): + mp_clump_slow(*args) + + context=multiprocessing.get_context('fork') + pool=context.Pool(processes=nproc, maxtasksperchild=1) + args=[] + for i in range(n_z): + args.append([i,obj_tmp[i], cy[i], cx[i], fwhm[i], sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose ]) + try: + pool.map_async(_mp_clump_slow, args, chunksize=1 ).get(timeout=10_000_000) + finally: + pool.close() + pool.join() + bpix_map_cumul = np.zeros_like(obj_tmp, dtype=obj_tmp.dtype) + bpix_map_cumul[:]=bpix_map_cumul_shared[:] + obj_tmp[:]=obj_tmp_shared_clump[:] + shm_clump.close() + shm_clump.unlink() + shm_clump_bpix.close() + shm_clump_bpix.unlink() + else: if isinstance(fwhm, (float, int)): fwhm_round = int(round(fwhm)) @@ -824,15 +869,50 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, fwhm_round = fwhm_round+1-(fwhm_round % 2) # make it odd neighbor_box = max(3, fwhm_round) # to not replace a companion nneig = sum(np.arange(3, neighbor_box+2, 2)) - for i in range(n_z): - if verbose: - print('************Frame # ', i, ' *************') - if bpm_mask.ndim == 3: - bpm = bpm_mask[i] - else: - bpm = bpm_mask - obj_tmp[i] = sigma_filter(obj_tmp[i], bpm, neighbor_box, - nneig, half_res_y, verbose) + + if nproc==1: + for i in range(n_z): + if verbose: + print('Using serial approach') + print('************Frame # ', i, ' *************') + if bpm_mask.ndim == 3: + bpm = bpm_mask[i] + else: + bpm = bpm_mask + obj_tmp[i] = sigma_filter(obj_tmp[i], bpm, neighbor_box, + nneig, half_res_y, verbose) + else: + msg="Cleaning frames using ADACS' multiprocessing appraoch" + print(msg) + #creating shared memory that each process writes into. + shm_clump = shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) + #creating an array that uses shared memory buffer and has the properties of obj_tmp. + obj_tmp_shared_clump = np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump.buf) + #function that is called repeatedly by each process. + def mp_clean_clump(j, obj_tmp, bpm, neighbor_box, nneig, half_res_y, verbose): + obj_tmp_shared_clump[j] = sigma_filter(obj_tmp, bpm, neighbor_box, nneig, half_res_y, verbose) + + global _mp_clean_clump + #function that converts the args into bite-sized pieces for mp_clean_clump. + def _mp_clean_clump(args): + mp_clean_clump(*args) + context=multiprocessing.get_context('fork') + pool=context.Pool(processes=nproc, maxtasksperchild=1) + args=[] + for j in range(n_z): + if bpm_mask.ndim == 3: + bpm = bpm_mask[j] + else: + bpm = bpm_mask + args.append([j,obj_tmp[j], bpm, neighbor_box, nneig, half_res_y, verbose]) + try: + pool.map_async(_mp_clean_clump, args, chunksize=1 ).get(timeout=10_000_000) + finally: + pool.close() + pool.join() + obj_tmp[:]=obj_tmp_shared_clump[:] + shm_clump.close() + shm_clump.unlink() bpix_map_cumul = bpm_mask # make it a binary map From 325a3ef737cbfe696f67f838f95a1b41fbc92596 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Wed, 14 Dec 2022 15:55:05 +0800 Subject: [PATCH 03/12] cube_fix_badpix_isolated function updated with multiprocessing --- vip_hci/preproc/badpixremoval.py | 80 ++++++++++++++++++++++++++------ 1 file changed, 65 insertions(+), 15 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index bd70bfc8..8c21c3d4 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -163,7 +163,7 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, sigma_clip=3, num_neig=5, size=5, frame_by_frame=False, protect_mask=0, cxy=None, mad=False, ignore_nan=True, verbose=True, - full_output=False): + full_output=False, nproc=1): """ Corrects the bad pixels, marked in the bad pixel mask. The bad pixel is replaced by the median of the adjacent pixels. This function is very fast but works only with isolated (sparse) pixels. @@ -213,6 +213,10 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, full_output: bool, {False,True}, optional Whether to return as well the cube of bad pixel maps and the cube of defined annuli. + nproc: int, optional + This feature is added following ADACS update. Refers to the number of processors + available for calculations. Choosing a number >1 enables multiprocessing for the + correction of frames. This happens only when frame_by_frame=True. Return ------ @@ -268,20 +272,66 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, if bpm_mask.ndim == 2: bpm_mask = [bpm_mask]*n_frames bpm_mask = np.array(bpm_mask) - for i in Progressbar(range(n_frames), desc="processing frames"): - if bpm_mask is not None: - bpm_mask_tmp = bpm_mask[i] - else: - bpm_mask_tmp = None - res = frame_fix_badpix_isolated(array[i], bpm_mask=bpm_mask_tmp, - sigma_clip=sigma_clip, - num_neig=num_neig, size=size, - protect_mask=protect_mask, - verbose=False, cxy=(cx[i], cy[i]), - ignore_nan=ignore_nan, - full_output=True) - array_out[i] = res[0] - final_bpm[i] = res[1] + if nproc==1: + for i in Progressbar(range(n_frames), desc="processing frames"): + if bpm_mask is not None: + bpm_mask_tmp = bpm_mask[i] + else: + bpm_mask_tmp = None + res = frame_fix_badpix_isolated(array[i], bpm_mask=bpm_mask_tmp, + sigma_clip=sigma_clip, + num_neig=num_neig, size=size, + protect_mask=protect_mask, + verbose=False, cxy=(cx[i], cy[i]), + ignore_nan=ignore_nan, + full_output=True) + array_out[i] = res[0] + final_bpm[i] = res[1] + else: + print("Processing using ADACS' multiprocessing approach...") + # multiprocessing included only in the frame-by-frame branch of the if statement above. + # creating shared memory buffer for the cube (array) + shm_array_out = shared_memory.SharedMemory(create=True, size=array.nbytes) + # creating a shared array_out version that is the shm_array_out buffer above. + shared_array_out = np.ndarray(array.shape, dtype=array.dtype, buffer=shm_array_out.buf) + # creating shared memory buffer for the final bad pixel mask cube. + shm_final_bpm = shared_memory.SharedMemory(create=True, size=final_bpm.nbytes) + # creating a shared final_bpm version that is in the shm_final_bpm buffer above. + shared_final_bpm = np.ndarray(final_bpm.shape, dtype=final_bpm.dtype, buffer=shm_final_bpm.buf) + + #function that calls frame_fix_badpix_isolated using the similar arguments as in if nproc==1 branch above. + def mp_clean_isolated (j,frame, bpm_mask=None, sigma_clip=3, num_neig=5, size=5, protect_mask=0, verbose=False, cxy=None, ignore_nan=True, full_output=True): + shared_array_out[j], shared_final_bpm[j] = frame_fix_badpix_isolated(frame, bpm_mask=bpm_mask, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, verbose=verbose, cxy=cxy, ignore_nan=ignore_nan, full_output=full_output) + #function that unwraps the arguments and passes them to mp_clean_isolated. + global _mp_clean_isolated + def _mp_clean_isolated (args): + pargs=args[0:2] + kwargs=args[2] + mp_clean_isolated(*pargs, **kwargs) + + context=multiprocessing.get_context('fork') + pool=context.Pool(processes=nproc, maxtasksperchild=1) + + args=[] + for j in range(n_frames): + if bpm_mask is not None: + bpm_mask_tmp = bpm_mask[j] + else: + bpm_mask_tmp = None + dict_kwargs={'bpm_mask' : bpm_mask_tmp, 'sigma_clip': sigma_clip, 'num_neig': num_neig, 'size' : size, 'protect_mask': protect_mask, 'cxy' : (cx[j], cy[j]), 'ignore_nan': ignore_nan} + args.append([j, array[j], dict_kwargs]) + + try: + pool.map_async(_mp_clean_isolated, args, chunksize=1 ).get(timeout=10_000_000) + finally: + pool.close() + pool.join() + array_out[:]=shared_array_out[:] + final_bpm[:]=shared_final_bpm[:] + shm_array_out.close() + shm_array_out.unlink() + shm_final_bpm.close() + shm_final_bpm.unlink() count_bp = np.sum(final_bpm) else: if bpm_mask is None or not correct_only: From e7102b982802408f5daa724aacdf728f06a6e843 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Tue, 3 Jan 2023 10:54:59 +0800 Subject: [PATCH 04/12] Using Iains updated clip_sigma --- vip_hci/stats/clip_sigma.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/vip_hci/stats/clip_sigma.py b/vip_hci/stats/clip_sigma.py index b9744f9f..59f2dd9a 100644 --- a/vip_hci/stats/clip_sigma.py +++ b/vip_hci/stats/clip_sigma.py @@ -287,7 +287,7 @@ def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, if bpm_mask_ori is None: gpm_ori = np.ones(array.shape) else: - gpm_ori = np.ones(array.shape)-bpm_mask_ori + gpm_ori = np.ones(array.shape) - bpm_mask_ori ny, nx = array.shape bpm = np.ones(array.shape) @@ -331,7 +331,11 @@ def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, x-hbox_l:x+hbox_r+1] gp_arr = gpm_ori[y-hbox_b:y+hbox_t+1, x-hbox_l:x+hbox_r+1] - neighbours = sub_arr[np.where(gp_arr)].flatten() + gp_idx = np.nonzero(gp_arr) + neighbours = [] + for n, (i, j) in enumerate(zip(gp_idx[0], gp_idx[1])): + neighbours.append(sub_arr[i, j]) + neighbours = np.array(neighbours) neigh_list = [] remove_itself = True for i in range(neighbours.shape[0]): From 350a1c96a959f813ce83f36dbdad80a1b3005f99 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Wed, 4 Jan 2023 11:33:57 +0800 Subject: [PATCH 05/12] Dummy call added on cube_fix_badpix_isolated and cube_correct_nan functions --- vip_hci/preproc/badpixremoval.py | 8 +++++++- vip_hci/preproc/cosmetics.py | 3 +++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 8c21c3d4..85520ae6 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -289,6 +289,13 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, final_bpm[i] = res[1] else: print("Processing using ADACS' multiprocessing approach...") + #dummy calling the function to create cached version of the code prior to forking + if bpm_mask is not None: + bpm_mask_dum = bpm_mask[0] + else: + bpm_mask_dum = None + #point of dummy call + dummy_call = frame_fix_badpix_isolated(array[0], bpm_mask=bpm_mask_dum, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, verbose=False, cxy=(cx[0], cy[0]), ignore_nan=ignore_nan, full_output=False) # multiprocessing included only in the frame-by-frame branch of the if statement above. # creating shared memory buffer for the cube (array) shm_array_out = shared_memory.SharedMemory(create=True, size=array.nbytes) @@ -884,7 +891,6 @@ def bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, shm_clump_bpix= shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) #works with dtype=obj_tmp.dtype but not dtype=int bpix_map_cumul_shared= np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump_bpix.buf) - def mp_clump_slow(j, obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose): obj_tmp_shared_clump[j], bpix_map_cumul_shared[j] = bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose) diff --git a/vip_hci/preproc/cosmetics.py b/vip_hci/preproc/cosmetics.py index f454bf12..4bbb0d06 100644 --- a/vip_hci/preproc/cosmetics.py +++ b/vip_hci/preproc/cosmetics.py @@ -348,6 +348,9 @@ def cube_correct_nan(cube, neighbor_box=3, min_neighbors=3, verbose=False, msg = "In channel {}, {} NaN pixels were corrected" print(msg.format(zz, nnanpix)) else: + #dummy calling the function to prevent compiling of the function by individual workers + #This should save some time. + dummy_obj = nan_corr_2d(obj_tmp[0], neighbor_box, min_neighbors, half_res_y, verbose, full_output=False) if verbose: msg = "Correcting NaNs in multiprocessing using ADACS' approach..." print(msg) From 1a25aa8d1e95cdff35c537a9bfb5620210ed95f2 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Mon, 9 Jan 2023 10:48:07 +0800 Subject: [PATCH 06/12] Dummy clump added on clump function for the branch correct_only=True --- vip_hci/preproc/badpixremoval.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 85520ae6..f78f12f3 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -940,6 +940,13 @@ def _mp_clump_slow(args): else: msg="Cleaning frames using ADACS' multiprocessing appraoch" print(msg) + #dummy calling sigma_filter function to create a cached version of the numba function + if bpm_mask.ndim == 3: + dummy_bpm = bpm_mask[i] + else: + dummy_bpm = bpm_mask + #Actual dummy call is here. + dummy_obj_tmp = sigma_filter(obj_tmp[0], dummy_bpm, neighbor_box, nneig, half_res_y, verbose) #creating shared memory that each process writes into. shm_clump = shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) #creating an array that uses shared memory buffer and has the properties of obj_tmp. From d9966af4657d41b33d692d42ef55fc0b8bea3306 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Fri, 13 Jan 2023 12:27:13 +0800 Subject: [PATCH 07/12] Badframes edited according to updated vip as of 13Jan2023 --- vip_hci/preproc/badframes.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/vip_hci/preproc/badframes.py b/vip_hci/preproc/badframes.py index c54285ac..633e660d 100644 --- a/vip_hci/preproc/badframes.py +++ b/vip_hci/preproc/badframes.py @@ -63,9 +63,14 @@ def cube_detect_badfr_pxstats(array, mode='annulus', in_radius=10, width=10, """ check_array(array, 3, msg='array') - if in_radius+width > array[0].shape[0]/2: - msgve = 'Inner radius and annulus size are too big (out of boundaries)' - raise ValueError(msgve) + if mode == 'annulus': + if in_radius + width > array[0].shape[0] / 2: + msgve = 'Inner radius and annulus size are too big (out of boundaries)' + raise ValueError(msgve) + elif mode == 'circle': + if in_radius > array[0].shape[0] / 2: + msgve = 'Radius size is too big (out of boundaries)' + raise ValueError(msgve) if verbose: start_time = time_ini() @@ -367,4 +372,4 @@ def cube_detect_badfr_correlation(array, frame_ref, crop_size=30, if full_output: return good_index_list, bad_index_list, distances else: - return good_index_list, bad_index_list + return good_index_list, bad_index_list \ No newline at end of file From e5551ec3daa8b82a120fe854b59f69d7fd68d640 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Fri, 13 Jan 2023 14:34:25 +0800 Subject: [PATCH 08/12] Detection.py updated from VIP master branch --- vip_hci/metrics/detection.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/vip_hci/metrics/detection.py b/vip_hci/metrics/detection.py index dc2813b7..1e56e2c8 100644 --- a/vip_hci/metrics/detection.py +++ b/vip_hci/metrics/detection.py @@ -322,6 +322,7 @@ def print_abort(): yy_out = [] xx_out = [] snr_list = [] + snr_final = [] if mode in ('lpeaks', 'log', 'dog'): xx -= pad @@ -342,6 +343,7 @@ def print_abort(): _ = frame_report(array, fwhm, (x, y), verbose=verbose) yy_final.append(y) xx_final.append(x) + snr_final.append(snr_value) else: yy_out.append(y) xx_out.append(x) @@ -353,17 +355,18 @@ def print_abort(): if verbose: print(sep) - if debug or full_output: + if debug: table_full = pn.DataFrame({'y': yy.tolist(), 'x': xx.tolist(), 'px_snr': snr_list}) table_full.sort_values('px_snr') + print(table_full) yy_final = np.array(yy_final) xx_final = np.array(xx_final) yy_out = np.array(yy_out) xx_out = np.array(xx_out) - table = pn.DataFrame({'y': yy_final.tolist(), 'x': xx_final.tolist()}) + table = pn.DataFrame({'y': yy_final.tolist(), 'x': xx_final.tolist(), 'px_snr': snr_final}) if plot: coords = tuple(zip(xx_out.tolist() + xx_final.tolist(), @@ -373,9 +376,6 @@ def print_abort(): plot_frames(array, dpi=120, circle=coords, circle_alpha=circlealpha, circle_label=True, circle_radius=fwhm, **kwargs) - if debug: - print(table_full) - if full_output: return table else: @@ -561,4 +561,4 @@ def mask_sources(mask, ap_rad): rad_arr = dist_matrix(ny, cx=s_coords[1][s], cy=s_coords[0][s]) mask_out[np.where(rad_arr < ap_rad)] = 0 - return mask_out + return mask_out \ No newline at end of file From 02c725f7f48243ddedbe7f6278512632b2ff7c13 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Fri, 13 Jan 2023 14:50:09 +0800 Subject: [PATCH 09/12] snr_source updated from VIP master --- vip_hci/metrics/snr_source.py | 99 +++++++++++++++++++++++++---------- 1 file changed, 70 insertions(+), 29 deletions(-) diff --git a/vip_hci/metrics/snr_source.py b/vip_hci/metrics/snr_source.py index f344d8d0..a7b65f12 100644 --- a/vip_hci/metrics/snr_source.py +++ b/vip_hci/metrics/snr_source.py @@ -9,6 +9,7 @@ __author__ = 'Carlos Alberto Gomez Gonzalez, O. Absil @ ULg, V. Christiaens' __all__ = ['snr', 'snrmap', + 'indep_ap_centers', 'significance', 'frame_report'] @@ -216,9 +217,69 @@ def _snr_approx(array, source_xy, fwhm, centery, centerx): snr_value = signal / noise return sourcey, sourcex, snr_value +def indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes=False, + exclude_theta_range=None): + + sourcex, sourcey = source_xy + centery, centerx = frame_center(array) + sep = dist(centery, centerx, float(sourcey), float(sourcex)) + theta_0 = np.rad2deg(np.arctan2(sourcey - centery, sourcex - centerx)) + + if exclude_theta_range is not None: + exc_theta_range = list(exclude_theta_range) + + if not sep > (fwhm / 2) + 1: + raise RuntimeError('`source_xy` is too close to the frame center') + + # sens = 'clock' # counterclock + # assumes clockwise rotation when building test apertures + # change sign and conditions if counterclockwise + sign = -1 + if exclude_theta_range is not None: + if theta_0 > exc_theta_range[0] and theta_0 < exc_theta_range[1]: + exc_theta_range[0] += 360 + while theta_0 < exc_theta_range[1]: + theta_0 += 360 + theta = theta_0 + + angle = np.arcsin(fwhm / 2. / sep) * 2 + number_apertures = int(np.floor(2 * np.pi / angle)) + yy = [] + xx = [] + yy_all = np.zeros(number_apertures) + xx_all = np.zeros(number_apertures) + cosangle = np.cos(angle) + sinangle = np.sin(angle) + xx.append(sourcex - centerx) + yy.append(sourcey - centery) + xx_all[0] = sourcex - centerx + yy_all[0] = sourcey - centery + + for i in range(number_apertures - 1): + xx_all[i + 1] = cosangle * xx_all[i] - sign * sinangle * yy_all[i] + yy_all[i + 1] = cosangle * yy_all[i] + sign * sinangle * xx_all[i] + theta += sign * np.rad2deg(angle) + if exclude_negative_lobes and (i == 0 or i == number_apertures - 2): + continue + if exclude_theta_range is None: + xx.append(cosangle * xx_all[i] - sign * sinangle * yy_all[i]) + yy.append(cosangle * yy_all[i] + sign * sinangle * xx_all[i]) + else: + if theta < exc_theta_range[0] or theta > exc_theta_range[1]: + xx.append(cosangle * xx_all[i] - sign * sinangle * yy_all[i]) + yy.append(cosangle * yy_all[i] + sign * sinangle * xx_all[i]) + + xx = np.array(xx) + yy = np.array(yy) + + xx += centerx + yy += centery + + return yy, xx def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False, - exclude_negative_lobes=False, plot=False, verbose=False): + exclude_negative_lobes=False, exclude_theta_range=None, plot=False, + verbose=False): """ Calculate the S/N (signal to noise ratio) of a test resolution element in a residual frame (e.g. post-processed with LOCI, PCA, etc). Implements @@ -261,6 +322,11 @@ def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False, exclude_negative_lobes : bool, opt Whether to include the adjacent aperture lobes to the tested location or not. Can be set to True if the image shows significant neg lobes. + exclude_theta_range : tuple of 2 floats or None, opt + If provided, range of trigonometric angles in deg (measured from + positive x axis), to be avoided for apertures used for noise estimation. + WARNING: this is to be used wisely, e.g. only if a known authentic + circumstellar signal is biasing the SNR estimate. plot : bool, optional Plots the frame and the apertures considered for clarity. verbose: bool, optional @@ -289,36 +355,11 @@ def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False, raise TypeError('`array2` has not the same shape as input array') sourcex, sourcey = source_xy - centery, centerx = frame_center(array) - sep = dist(centery, centerx, float(sourcey), float(sourcex)) - - if not sep > (fwhm/2)+1: - raise RuntimeError('`source_xy` is too close to the frame center') - - sens = 'clock' # counterclock - angle = np.arcsin(fwhm/2./sep)*2 - number_apertures = int(np.floor(2*np.pi/angle)) - yy = np.zeros((number_apertures)) - xx = np.zeros((number_apertures)) - cosangle = np.cos(angle) - sinangle = np.sin(angle) - xx[0] = sourcex - centerx - yy[0] = sourcey - centery - for i in range(number_apertures-1): - if sens == 'clock': - xx[i+1] = cosangle*xx[i] + sinangle*yy[i] - yy[i+1] = cosangle*yy[i] - sinangle*xx[i] - elif sens == 'counterclock': - xx[i+1] = cosangle*xx[i] - sinangle*yy[i] - yy[i+1] = cosangle*yy[i] + sinangle*xx[i] + yy, xx = indep_ap_centers(array, source_xy, fwhm, exclude_negative_lobes, + exclude_theta_range) - xx += centerx - yy += centery rad = fwhm/2. - if exclude_negative_lobes: - xx = np.concatenate(([xx[0]], xx[2:-1])) - yy = np.concatenate(([yy[0]], yy[2:-1])) apertures = photutils.CircularAperture( zip(xx, yy), r=rad) # Coordinates (X,Y) @@ -532,4 +573,4 @@ def frame_report(array, fwhm, source_xy=None, verbose=True): print(msg3.format(stdsnr_i)) print(sep) - return source_xy, obj_flux, snr_centpx, meansnr_pixels + return source_xy, obj_flux, snr_centpx, meansnr_pixels \ No newline at end of file From 6d4cbf792e75210b7ad3cdc9b9c68cf2b45a54f4 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Sat, 14 Jan 2023 09:28:55 +0800 Subject: [PATCH 10/12] recentering and rescaling patched from main VIP branch --- vip_hci/preproc/recentering.py | 5 +++-- vip_hci/preproc/rescaling.py | 22 ++++++++++++++++++---- 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/vip_hci/preproc/recentering.py b/vip_hci/preproc/recentering.py index 25da3466..214e85cd 100644 --- a/vip_hci/preproc/recentering.py +++ b/vip_hci/preproc/recentering.py @@ -1476,7 +1476,8 @@ def cube_recenter_2dfit(array, xy=None, fwhm=4, subi_size=5, model='gauss', if nproc == 1: res = [] - print('2d {}-fitting'.format(model)) + if verbose: + print('2d {}-fitting'.format(model)) for i in Progressbar(range(n_frames), desc="frames", verbose=verbose): if model == "2gauss": args = [array, i, subi_size, pos_y, pos_x, debug, fwhm[i], @@ -2027,4 +2028,4 @@ def _fit_2dannulus(array, fwhm=4, crop=False, cent=None, cropsize=15, return mean_y, mean_x else: final_hole_rad = best_rad[i_max, j_max]/fwhm - return mean_y, mean_x, final_hole_rad + return mean_y, mean_x, final_hole_rad \ No newline at end of file diff --git a/vip_hci/preproc/rescaling.py b/vip_hci/preproc/rescaling.py index e19a5280..d213805d 100644 --- a/vip_hci/preproc/rescaling.py +++ b/vip_hci/preproc/rescaling.py @@ -23,7 +23,7 @@ from scipy.ndimage import geometric_transform, zoom from scipy.optimize import minimize -from ..var import frame_center, get_square +from ..var import frame_center, get_square, cube_filter_highpass from .subsampling import cube_collapse from .recentering import frame_shift from .cosmetics import frame_crop @@ -693,7 +693,7 @@ def check_scal_vector(scal_vec): def find_scal_vector(cube, lbdas, fluxes, mask=None, nfp=2, fm="stddev", simplex_options=None, debug=False, imlib='vip-fft', - interpolation='lanczos4', **kwargs): + interpolation='lanczos4', hpf=False, fwhm_max=5, **kwargs): """ Find the optimal scaling factor for the channels of an IFS cube (or of dual-band pairs of images). @@ -723,6 +723,13 @@ def find_scal_vector(cube, lbdas, fluxes, mask=None, nfp=2, fm="stddev", pixels. options: dict, optional The scipy.optimize.minimize options. + hpf: bool, optional + Whether to high-pass filter images before searching for optimal scaling + factors. Can help to not be biased by a diffuse halo, and just leverage + speckle expansion. + max_fwhm: float, optional + Maximum FWHM of the PSF across all wavelengths, in pixels. Only used if + hpf is set to True. **kwargs: optional Optional arguments to the scipy.optimize.minimize function @@ -746,9 +753,16 @@ def find_scal_vector(cube, lbdas, fluxes, mask=None, nfp=2, fm="stddev", 'maxfev': 2000} scal_vec = np.ones(n_z) flux_vec = np.ones(n_z) + array = cube.copy() + if hpf: + med_sz = int(5*fwhm_max) + if not med_sz%2: + med_sz+=1 + array = cube_filter_highpass(cube, mode='median-subt', + median_size=med_sz) for z in range(n_z-1): flux_scal = fluxes[-1]/fluxes[z] - cube_tmp = np.array([cube[z], cube[-1]]) + cube_tmp = np.array([array[z], array[-1]]) if nfp == 1: p_ini = (scal_vec_ini[z],) solu = minimize(_chisquare_scal, p_ini, args=(cube_tmp, flux_scal, @@ -1050,4 +1064,4 @@ def scale_fft(array, scale, ori_dim=False): scaled[-kf_io:-kf_io+dim_pp, -kf_io:-kf_io+dim_pp] = array_resc array_resc = scaled - return array_resc + return array_resc \ No newline at end of file From a28fa44f464550c566f24f10595bd117ebac0d0c Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Mon, 16 Jan 2023 16:16:54 +0800 Subject: [PATCH 11/12] renamed variables after merge --- vip_hci/preproc/badpixremoval.py | 44 ++++++++++++++++---------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 4d7c4544..434d3c16 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -295,7 +295,7 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, correct_only=False, else: bpm_mask_dum = None #point of dummy call - dummy_call = frame_fix_badpix_isolated(array[0], bpm_mask=bpm_mask_dum, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, verbose=False, cxy=(cx[0], cy[0]), ignore_nan=ignore_nan, full_output=False) + frame_fix_badpix_isolated(array[0], bpm_mask=bpm_mask_dum, sigma_clip=sigma_clip, num_neig=num_neig, size=size, protect_mask=protect_mask, verbose=False, cxy=(cx[0], cy[0]), ignore_nan=ignore_nan, full_output=False) # multiprocessing included only in the frame-by-frame branch of the if statement above. # creating shared memory buffer for the cube (array) shm_array_out = shared_memory.SharedMemory(create=True, size=array.nbytes) @@ -869,11 +869,11 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, if isinstance(fwhm, (float, int)): fwhm = [fwhm]*n_z if nproc==1: - bpix_map_cumul = np.zeros_like(obj_tmp) + bpix_map_cumul = np.zeros_like(array_corr) for i in range(n_z): if verbose: print('************Frame # ', i, ' *************') - obj_tmp[i], bpix_map_cumul[i] = bp_removal_2d(obj_tmp[i], cy[i], + array_corr[i], bpix_map_cumul[i] = bp_removal_2d(array_corr[i], cy[i], cx[i], fwhm[i], sig, protect_mask, bpm_mask, min_thr, @@ -883,14 +883,14 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, msg="Cleaning frames using ADACS' multiprocessing appraoch" print(msg) #creating shared memory buffer space for the image cube. - shm_clump= shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) - obj_tmp_shared_clump = np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump.buf) + shm_clump= shared_memory.SharedMemory(create=True, size=array_corr.nbytes) + obj_tmp_shared_clump = np.ndarray(array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump.buf) #creating shared memory buffer space for the bad pixel cube. - shm_clump_bpix= shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) + shm_clump_bpix= shared_memory.SharedMemory(create=True, size=array_corr.nbytes) #works with dtype=obj_tmp.dtype but not dtype=int - bpix_map_cumul_shared= np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump_bpix.buf) - def mp_clump_slow(j, obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose): - obj_tmp_shared_clump[j], bpix_map_cumul_shared[j] = bp_removal_2d(obj_tmp, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose) + bpix_map_cumul_shared= np.ndarray(array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump_bpix.buf) + def mp_clump_slow(j, array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose): + obj_tmp_shared_clump[j], bpix_map_cumul_shared[j] = bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose) global _mp_clump_slow @@ -901,15 +901,15 @@ def _mp_clump_slow(args): pool=context.Pool(processes=nproc, maxtasksperchild=1) args=[] for i in range(n_z): - args.append([i,obj_tmp[i], cy[i], cx[i], fwhm[i], sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose ]) + args.append([i,array_corr[i], cy[i], cx[i], fwhm[i], sig, protect_mask, bpm_mask, min_thr, half_res_y, mad, verbose ]) try: pool.map_async(_mp_clump_slow, args, chunksize=1 ).get(timeout=10_000_000) finally: pool.close() pool.join() - bpix_map_cumul = np.zeros_like(obj_tmp, dtype=obj_tmp.dtype) + bpix_map_cumul = np.zeros_like(array_corr, dtype=array_corr.dtype) bpix_map_cumul[:]=bpix_map_cumul_shared[:] - obj_tmp[:]=obj_tmp_shared_clump[:] + array_corr[:]=obj_tmp_shared_clump[:] shm_clump.close() shm_clump.unlink() shm_clump_bpix.close() @@ -933,25 +933,25 @@ def _mp_clump_slow(args): bpm = bpm_mask[i] else: bpm = bpm_mask - obj_tmp[i] = sigma_filter(obj_tmp[i], bpm, neighbor_box, + array_corr[i] = sigma_filter(array_corr[i], bpm, neighbor_box, nneig, half_res_y, verbose) else: msg="Cleaning frames using ADACS' multiprocessing appraoch" print(msg) #dummy calling sigma_filter function to create a cached version of the numba function if bpm_mask.ndim == 3: - dummy_bpm = bpm_mask[i] + dummy_bpm = bpm_mask[0] else: dummy_bpm = bpm_mask #Actual dummy call is here. - dummy_obj_tmp = sigma_filter(obj_tmp[0], dummy_bpm, neighbor_box, nneig, half_res_y, verbose) + sigma_filter(array_corr[0], dummy_bpm, neighbor_box, nneig, half_res_y, verbose) #creating shared memory that each process writes into. - shm_clump = shared_memory.SharedMemory(create=True, size=obj_tmp.nbytes) - #creating an array that uses shared memory buffer and has the properties of obj_tmp. - obj_tmp_shared_clump = np.ndarray(obj_tmp.shape, dtype=obj_tmp.dtype, buffer=shm_clump.buf) + shm_clump = shared_memory.SharedMemory(create=True, size=array_corr.nbytes) + #creating an array that uses shared memory buffer and has the properties of array_corr. + obj_tmp_shared_clump = np.ndarray(array_corr.shape, dtype=array_corr.dtype, buffer=shm_clump.buf) #function that is called repeatedly by each process. - def mp_clean_clump(j, obj_tmp, bpm, neighbor_box, nneig, half_res_y, verbose): - obj_tmp_shared_clump[j] = sigma_filter(obj_tmp, bpm, neighbor_box, nneig, half_res_y, verbose) + def mp_clean_clump(j, array_corr, bpm, neighbor_box, nneig, half_res_y, verbose): + obj_tmp_shared_clump[j] = sigma_filter(array_corr, bpm, neighbor_box, nneig, half_res_y, verbose) global _mp_clean_clump #function that converts the args into bite-sized pieces for mp_clean_clump. @@ -965,13 +965,13 @@ def _mp_clean_clump(args): bpm = bpm_mask[j] else: bpm = bpm_mask - args.append([j,obj_tmp[j], bpm, neighbor_box, nneig, half_res_y, verbose]) + args.append([j,array_corr[j], bpm, neighbor_box, nneig, half_res_y, verbose]) try: pool.map_async(_mp_clean_clump, args, chunksize=1 ).get(timeout=10_000_000) finally: pool.close() pool.join() - obj_tmp[:]=obj_tmp_shared_clump[:] + array_corr[:]=obj_tmp_shared_clump[:] shm_clump.close() shm_clump.unlink() bpix_map_cumul = bpm_mask From 5da8fdc3b5ba6bc3a49859a35969665b1a76dbd2 Mon Sep 17 00:00:00 2001 From: SrikanthKom <19237130@student.curtin.edu.au> Date: Thu, 19 Jan 2023 11:48:54 +0800 Subject: [PATCH 12/12] shared_memory issue resolved for python3.7 --- vip_hci/preproc/badpixremoval.py | 13 ++++++++++++- vip_hci/preproc/cosmetics.py | 15 +++++++++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/vip_hci/preproc/badpixremoval.py b/vip_hci/preproc/badpixremoval.py index 434d3c16..5c67fe4b 100644 --- a/vip_hci/preproc/badpixremoval.py +++ b/vip_hci/preproc/badpixremoval.py @@ -36,7 +36,18 @@ from .cosmetics import frame_pad from multiprocessing import Process import multiprocessing -from multiprocessing import shared_memory, set_start_method +from multiprocessing import set_start_method +try: + from multiprocessing import shared_memory +except ImportError: + print('Failed to import shared_memory from multiprocessing') + try: + print('Trying to import shared_memory directly(for python 3.7)') + import shared_memory + except ModuleNotFoundError: + print('Use shared_memory on python 3.7 to activate') + print('multiprocessing on badpixels using..') + print('pip install shared-memory38') import warnings try: diff --git a/vip_hci/preproc/cosmetics.py b/vip_hci/preproc/cosmetics.py index 8c4eb4e6..7bd0d82c 100644 --- a/vip_hci/preproc/cosmetics.py +++ b/vip_hci/preproc/cosmetics.py @@ -1,5 +1,5 @@ #! /usr/bin/env python -# check123 + """ Module with cosmetics procedures. Contains the function for bad pixel fixing. Also functions for cropping cubes. @@ -23,7 +23,18 @@ from ..var import frame_center, get_square from multiprocessing import Process import multiprocessing -from multiprocessing import shared_memory, set_start_method +from multiprocessing import set_start_method +try: + from multiprocessing import shared_memory +except ImportError: + print('Failed to import shared_memory from multiprocessing') + try: + print('Trying to import shared_memory directly(for python 3.7)') + import shared_memory + except ModuleNotFoundError: + print('Use shared_memory on python 3.7 to activate') + print('multiprocessing on badpixels using..') + print('pip install shared-memory38') def cube_crop_frames(array, size, xy=None, force=False, verbose=True,