diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ccf4b1a4..92699db1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -38,18 +38,18 @@ jobs: run: | python -m pip install --upgrade pip pip install --editable . --group dev -# Uncomment below when ready to deal with a lot of PEP8 formatting changes -# - name: Verify files with pre-commit -# run: | -# # Setup pre-commit hooks -# pre-commit clean -# pre-commit install --hook-type pre-merge-commit -# pre-commit install --hook-type pre-push -# pre-commit install --hook-type post-rewrite -# pre-commit install-hooks -# pre-commit install -# # Run pre-commit hooks -# pre-commit run --files src/**/*.py + # Uncomment below when ready to deal with a lot of PEP8 formatting changes + - name: Verify files with pre-commit + run: | + # Setup pre-commit hooks + pre-commit clean + pre-commit install --hook-type pre-merge-commit + pre-commit install --hook-type pre-push + pre-commit install --hook-type post-rewrite + pre-commit install-hooks + pre-commit install + # Run pre-commit hooks + pre-commit run --files src/**/*.py - name: Test with pytest run: | coverage run --source src -m pytest tests/pre_3_10 --splits 3 --group ${{ matrix.group }} @@ -75,17 +75,17 @@ jobs: run: | python -m pip install --upgrade pip pip install --editable . --group dev -# - name: Verify files with pre-commit -# run: | -# # Setup pre-commit hooks -# pre-commit clean -# pre-commit install --hook-type pre-merge-commit -# pre-commit install --hook-type pre-push -# pre-commit install --hook-type post-rewrite -# pre-commit install-hooks -# pre-commit install -# # Run pre-commit hooks -# pre-commit run --files src/vip_hci/objects/*.py + - name: Verify files with pre-commit + run: | + # Setup pre-commit hooks + pre-commit clean + pre-commit install --hook-type pre-merge-commit + pre-commit install --hook-type pre-push + pre-commit install --hook-type post-rewrite + pre-commit install-hooks + pre-commit install + # Run pre-commit hooks + pre-commit run --files src/vip_hci/objects/*.py - name: Test with pytest run: | coverage run --source=src/vip_hci/objects/ -m pytest tests/post_3_10 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1071a2a0..91a09a0e 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,7 +6,7 @@ repos: - id: end-of-file-fixer - id: check-yaml - - repo: https://github.com/asottile/reorder_python_imports - rev: v3.9.0 - hooks: - - id: reorder-python-imports +# - repo: https://github.com/asottile/reorder_python_imports +# rev: v3.9.0 +# hooks: +# - id: reorder-python-imports diff --git a/docs/source/conf.py b/docs/source/conf.py index 7c16a33b..f2633656 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -53,6 +53,9 @@ suppress_warnings = ["mystnb.unknown_mime_type"] +# default timeout is 30s which is too short to compile most VIP notebooks +nb_execution_timeout = 99999 + # Add any paths that contain templates here, relative to this directory. templates_path = []#['./_templates'] diff --git a/docs/source/faq.rst b/docs/source/faq.rst index daca7397..54b9bcdf 100644 --- a/docs/source/faq.rst +++ b/docs/source/faq.rst @@ -96,4 +96,4 @@ your ``~/.bash_profile``: .. rubric:: Where is the ``specfit`` module? -From versions 1.0.0 to 1.0.3, ``specfit`` was a module of VIP offering atmosphere retrieval and spectral characterisation of directly imaged companions. Given the divergence with the original purpose of VIP, starting from version 1.1.0, it has been renamed, expanded, moved to a separate `GitHub repository `_ and converted into its own `package `_ (called ``special``). \ No newline at end of file +From versions 1.0.0 to 1.0.3, ``specfit`` was a module of VIP offering atmosphere retrieval and spectral characterisation of directly imaged companions. Given the divergence with the original purpose of VIP, starting from version 1.1.0, it has been renamed, expanded, moved to a separate `GitHub repository `_ and converted into its own `package `_ (called ``special``). diff --git a/docs/source/tutorials/02_preproc.ipynb b/docs/source/tutorials/02_preproc.ipynb index d42ba04b..a96ae756 100644 --- a/docs/source/tutorials/02_preproc.ipynb +++ b/docs/source/tutorials/02_preproc.ipynb @@ -1176,7 +1176,7 @@ "metadata": {}, "source": [ "```{note}\n", - "An alternative to the quick NaN replacement above, is to (i) create a binary bad pixel map which includes the NaN values (using e.g. `np.isnan`), (ii) add these bad pixels to the bad pixel map inferred through one or more methods detailed in the next section, and (iii) finally proceed with a single correction of all NaN and bad pixels together with the `cube_fix_badpix_interp` routine.\n", + "A better alternative to the quick NaN replacement above, is to (i) create a binary bad pixel map which includes the NaN values (using e.g. `np.isnan`), (ii) add these bad pixels to the bad pixel map inferred through one or more methods detailed in the next section, and (iii) finally proceed with a single correction of all NaN and bad pixels together with the `cube_fix_badpix_interp` routine.\n", "```" ] }, diff --git a/src/vip_hci/config/utils_param.py b/src/vip_hci/config/utils_param.py index 1f202f27..c7d9bc1c 100644 --- a/src/vip_hci/config/utils_param.py +++ b/src/vip_hci/config/utils_param.py @@ -8,7 +8,9 @@ KWARGS_EXCEPTIONS = ["param"] -def filter_duplicate_keys(filter_item: any, ref_item: any, filter_in: bool = True): +def filter_duplicate_keys(filter_item: any, + ref_item: any, + filter_in: bool = True): """ Filter in or out keys of an item based on a reference item. @@ -30,14 +32,16 @@ def filter_duplicate_keys(filter_item: any, ref_item: any, filter_in: bool = Tru elif isinstance(filter_item, object): filter_dict = vars(filter_item).copy() else: - raise TypeError("The item to be filtered is neither a dictionnary or an object") + msg = "The item to be filtered is neither a dictionnary or an object" + raise TypeError(msg) if isinstance(ref_item, dict): ref_dict = ref_item.copy() elif isinstance(ref_item, object): ref_dict = vars(ref_item).copy() else: - raise TypeError("The reference item is neither a dictionnary or an object") + msg = "The reference item is neither a dictionnary or an object" + raise TypeError(msg) # Find keys that must serve as a filter for `filter_item` common_keys = set(filter_dict.keys()) & set(ref_dict.keys()) @@ -64,9 +68,9 @@ def setup_parameters( """ Help creating a dictionnary of parameters for a given function. - Look for the exact list of parameters needed for the ``fkt`` function and takes - only the attributes needed from the ``params_obj``. More parameters can be - included with the ``**add_pararms`` dictionnary. + Look for the exact list of parameters needed for the ``fkt`` function and + take only the attributes needed from the ``params_obj``. More parameters can + be included with the ``**add_pararms`` dictionnary. Parameters ---------- @@ -83,7 +87,7 @@ def setup_parameters( Returns ------- params_setup : dictionnary or list - The dictionnary comprised of parameters needed for the function, selected + The dictionnary comprised of parameters needed for the function selected amongst attributes of PostProc objects and additionnal parameters. Can be a list if asked for (used in specific cases such as when calling functions through ``vip_hci.config.utils_conf.pool_map``, see an example @@ -93,7 +97,8 @@ def setup_parameters( wanted_params = OrderedDict(signature(fkt).parameters) # Remove dupe keys in params_obj from add_params if add_params is not None: - obj_params = filter_duplicate_keys(filter_item=params_obj, ref_item=add_params) + obj_params = filter_duplicate_keys(filter_item=params_obj, + ref_item=add_params) all_params = {**obj_params, **add_params} else: all_params = vars(params_obj) diff --git a/src/vip_hci/fm/negfc_fmerit.py b/src/vip_hci/fm/negfc_fmerit.py index a680057b..822a8e2d 100644 --- a/src/vip_hci/fm/negfc_fmerit.py +++ b/src/vip_hci/fm/negfc_fmerit.py @@ -816,47 +816,71 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm, **algo_opt_copy, ) - elif algo == pca_annular: + elif algo == pca_annular or algo == nmf_annular: tol = algo_opt_copy.pop("tol", 1e-1) min_frames_lib = algo_opt_copy.pop("min_frames_lib", 2) max_frames_lib = algo_opt_copy.pop("max_frames_lib", 200) - nproc = algo_opt_copy.pop("nproc", 1) + radius_int = max(1, int(np.floor(r_guess - annulus_width / 2))) + radius_int = algo_opt_copy.pop("radius_int", radius_int) + asize = algo_opt_copy.pop("asize", annulus_width) + delta_rot = algo_opt_copy.pop("delta_rot", delta_rot) + _ = algo_opt_copy.pop("verbose", verbose) # crop cube to just be larger than annulus => FASTER PCA - crop_sz = int(2 * np.ceil(radius_int + annulus_width + 1)) + crop_sz = int(2 * np.ceil(radius_int + asize + 1)) if not crop_sz % 2: crop_sz += 1 - if crop_sz < array.shape[1] and crop_sz < array.shape[2]: - pad = int((array.shape[1] - crop_sz) / 2) - crop_cube = cube_crop_frames(array, crop_sz, verbose=False) + if crop_sz < cube.shape[-2] and crop_sz < cube.shape[-1]: + pad = int((cube.shape[-2] - crop_sz) / 2) + crop_cube = cube_crop_frames(cube, crop_sz, verbose=False) else: + crop_cube = cube pad = 0 - crop_cube = array - - pca_res_tmp = pca_annular( - cube=crop_cube, - angle_list=angs, - radius_int=radius_int, - fwhm=fwhm, - asize=annulus_width, - delta_rot=delta_rot, - ncomp=ncomp, - svd_mode=svd_mode, - scaling=scaling, - imlib=imlib, - interpolation=interpolation, - collapse=collapse, - tol=tol, - nproc=nproc, - min_frames_lib=min_frames_lib, - max_frames_lib=max_frames_lib, - full_output=False, - verbose=False, - weights=weights, - **algo_opt_copy, - ) + if algo == pca_annular: + res_tmp = algo( + cube=crop_cube, + angle_list=angs, + cube_ref=cube_ref, + radius_int=radius_int, + fwhm=fwhm, + asize=asize, + delta_rot=delta_rot, + ncomp=ncomp, + svd_mode=svd_mode, + scaling=scaling, + imlib=imlib, + interpolation=interpolation, + collapse=collapse, + weights=weights, + tol=tol, + min_frames_lib=min_frames_lib, + max_frames_lib=max_frames_lib, + full_output=False, + verbose=False, + **algo_opt_copy, + ) + else: + res_tmp = algo( + cube=crop_cube, + angle_list=angs, + cube_ref=cube_ref, + radius_int=radius_int, + fwhm=fwhm, + asize=annulus_width, + delta_rot=delta_rot, + ncomp=ncomp, + scaling=scaling, + imlib=imlib, + interpolation=interpolation, + collapse=collapse, + weights=weights, + min_frames_lib=min_frames_lib, + max_frames_lib=max_frames_lib, + full_output=False, + verbose=False, + **algo_opt_copy, + ) # pad again now - pca_res = np.pad(pca_res_tmp, pad, mode="constant", - constant_values=0) + pca_res = np.pad(res_tmp, pad, mode="constant", constant_values=0) if f_guess is not None and psfn is not None: pca_res_tinv = pca_annular( @@ -873,7 +897,6 @@ def get_mu_and_sigma(cube, angs, ncomp, annulus_width, aperture_radius, fwhm, interpolation=interpolation, collapse=collapse, tol=tol, - nproc=nproc, min_frames_lib=min_frames_lib, max_frames_lib=max_frames_lib, full_output=False, diff --git a/src/vip_hci/invprob/__init__.py b/src/vip_hci/invprob/__init__.py index a6a7c7bd..7f0915f0 100644 --- a/src/vip_hci/invprob/__init__.py +++ b/src/vip_hci/invprob/__init__.py @@ -1,6 +1,6 @@ """ Subpackage ``invprob`` aims to contain post-processing algorithms based on an -inverse problem approach, such as ANDROMEDA [MUG09]_ / [CAN15]_, Foward Model +inverse problem approach, such as ANDROMEDA [MUG09]_ / [CAN15]_, Forward Model Matched Filter [RUF17]_ / [DAH21a]_ or PACO [FLA18]_. """ diff --git a/src/vip_hci/invprob/utils_andro.py b/src/vip_hci/invprob/utils_andro.py index d3d1c81d..0802d932 100644 --- a/src/vip_hci/invprob/utils_andro.py +++ b/src/vip_hci/invprob/utils_andro.py @@ -331,4 +331,4 @@ def subpixel_shift(image, xshift, yshift): image_ft = np.fft.fft2(image) # no np.fft.fftshift applied! shifted_image = np.fft.ifft2(image_ft * fact).real - return shifted_image \ No newline at end of file + return shifted_image diff --git a/src/vip_hci/metrics/completeness.py b/src/vip_hci/metrics/completeness.py index 4a481321..27b8a515 100644 --- a/src/vip_hci/metrics/completeness.py +++ b/src/vip_hci/metrics/completeness.py @@ -166,7 +166,7 @@ def _estimate_snr_fc( mask = get_annulus_segments(frame_fin, (fwhm_med / 2) + 2, width, mode="mask")[ 0 ] - bmask = np.ma.make_mask(mask) + bmask = np.ma.make_mask(mask, shrink=False) yy, xx = np.where(bmask) if approximated: diff --git a/src/vip_hci/metrics/snr_source.py b/src/vip_hci/metrics/snr_source.py index ba277027..b651c0dd 100644 --- a/src/vip_hci/metrics/snr_source.py +++ b/src/vip_hci/metrics/snr_source.py @@ -85,7 +85,7 @@ def snrmap(array, fwhm, approximated=False, plot=False, known_sources=None, snrmap_array = np.zeros_like(array) width = min(sizey, sizex) / 2 - 1.5*fwhm mask = get_annulus_segments(array, fwhm, width, mode="mask")[0] - mask = np.ma.make_mask(mask) + mask = np.ma.make_mask(mask, shrink=False) # by making a bool mask *after* applying the mask to the array, we also mask # out zero values from the array. This logic cannot be simplified by using # mode="ind"! @@ -105,25 +105,28 @@ def snrmap(array, fwhm, approximated=False, plot=False, known_sources=None, width = min(sizey, sizex) / 2 - 1.5 * fwhm mask = get_annulus_segments(array, (fwhm / 2) + 1, width - 1, mode="mask")[0] - mask = np.ma.make_mask(mask) + mask = np.ma.make_mask(mask, shrink=False) yy, xx = np.where(mask) coords = [(int(x), int(y)) for (x, y) in zip(xx, yy)] res = pool_map(nproc, _snr_approx, array, iterable(coords), fwhm, cy, cx) - res = np.array(res, dtype=object) - yy = res[:, 0] - xx = res[:, 1] - snr_value = res[:, 2] + #res = np.array(res, dtype=object) + yy = np.array([res[i][0] for i in range(len(res))]) + xx = np.array([res[i][1] for i in range(len(res))]) + snr_value = np.array([res[i][2] for i in range(len(res))]) snrmap_array[yy.astype(int), xx.astype(int)] = snr_value # computing s/n map with Mawet+14 definition else: res = pool_map(nproc, snr, array, iterable(coords), fwhm, True, array2, use2alone, exclude_negative_lobes) - res = np.array(res, dtype=object) - yy = res[:, 0] - xx = res[:, 1] - snr_value = res[:, -1] + #res = np.array(res, dtype=object) + yy = np.array([res[i][0] for i in range(len(res))]) + xx = np.array([res[i][1] for i in range(len(res))]) + snr_value = np.array([res[i][-1] for i in range(len(res))]) + #yy = res[:][0] + #xx = res[:][1] + #snr_value = res[:][-1] snrmap_array[yy.astype('int'), xx.astype('int')] = snr_value # masking known sources diff --git a/src/vip_hci/preproc/badpixremoval.py b/src/vip_hci/preproc/badpixremoval.py index 5894a455..cecc3989 100644 --- a/src/vip_hci/preproc/badpixremoval.py +++ b/src/vip_hci/preproc/badpixremoval.py @@ -837,14 +837,15 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, cx=None, fwhm=4., sig=4., protect_mask=0, excl_mask=None, half_res_y=False, min_thr=None, max_nit=15, mad=True, bad_values=None, verbose=True, - full_output=False, nproc=1): + full_output=False, debug=True, 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 - clumps will be searched iteratively and replaced by the median of good - neighbouring pixel values, when enough of them are available. The size of - the box is set by the closest odd integer larger than fwhm (to avoid - accidentally replacing point sources). + 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 clumps will be searched iteratively and replaced by + the median of good neighbouring pixel values, when enough of them are + available. The size of the box is set by the closest odd integer larger than + fwhm (to avoid accidentally replacing point sources). Parameters ---------- @@ -910,6 +911,9 @@ 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. + debug: bool, {False,True}, optional + In case the algorithm encounters a problem for a given frame of the + cube, whether to stop the function and write a fits file of it. nproc: int, optional This feature is added following ADACS update. Refers to the number of processors available for calculations. Choosing a number >1 enables @@ -921,8 +925,8 @@ def cube_fix_badpix_clump(array, bpm_mask=None, correct_only=False, cy=None, The bad pixel corrected frame/cube. bpix_map: 2d or 3d array [full_output=True] The bad pixel map or the cube of bpix maps - """ + """ array_corr = array.copy() ndims = array_corr.ndim assert ndims == 2 or ndims == 3, "Object is not two or three dimensional.\n" @@ -995,9 +999,17 @@ def bp_removal_2d(array_corr, cy, cx, fwhm, sig, protect_mask, bpm_mask_ori, bpm_mask = excl_mask.copy() if bpm_mask_ori is not None: bpm_mask += bpm_mask_ori.astype(bool) - bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, - neighbor=True, num_neighbor=neighbor_box, mad=mad, - half_res_y=half_res_y) + try: + bp = clip_array(array_corr, sig, sig, bpm_mask, out_good=False, + neighbor=True, num_neighbor=neighbor_box, mad=mad, + half_res_y=half_res_y) + except AssertionError: + msg = "Prob with clip array using lower and upper sigma set to {}," + msg += "{} bad pixels in in input bad pixel mask," + msg += "{} neighbours in neighbour box" + print(msg.format(sig, np.sum(bpm_mask), neighbor_box)) + from vip_hci.fits import write_fits + write_fits("TMP.fits", array_corr) bpix_map = np.zeros_like(array_corr) bpix_map[bp] = 1 if min_thr is not None: diff --git a/src/vip_hci/preproc/recentering.py b/src/vip_hci/preproc/recentering.py index 0ef5e99d..5b7fff73 100644 --- a/src/vip_hci/preproc/recentering.py +++ b/src/vip_hci/preproc/recentering.py @@ -360,6 +360,10 @@ def frame_center_satspots(array, xy, subi_size=19, sigfactor=6, shift=False, beyond the edge with zeros. With 'mirror', the input is extended by reflecting about the center of the last pixel. With 'wrap', the input is extended by wrapping around to the opposite edge. Default is 'reflect'. + imlib : str, optional + See the documentation of the ``vip_hci.preproc.frame_shift`` function. + interpolation : str, optional + See the documentation of the ``vip_hci.preproc.frame_shift`` function. debug : bool, optional If True debug information is printed and plotted. verbose : bool, optional @@ -517,9 +521,10 @@ def intersection(L1, L2): def cube_recenter_satspots(array, xy, subi_size=19, sigfactor=6, plot=True, fit_type='moff', lbda=None, filter_freq=(0, 0), - border_mode='constant', debug=False, verbose=True, + border_mode='constant', imlib='vip-fft', + interpolation='lanczos4', debug=False, verbose=True, full_output=False): - """Recenter an image cube based on satellite spots (more details in `. + """Recenter an image cube based on satellite spots. The function relies on ``frame_center_satspots`` to align each image of the cube individually (details in ``vip_hci.preproc.frame_center_satspots``). @@ -572,6 +577,10 @@ def cube_recenter_satspots(array, xy, subi_size=19, sigfactor=6, plot=True, beyond the edge with zeros. With 'mirror', the input is extended by reflecting about the center of the last pixel. With 'wrap', the input is extended by wrapping around to the opposite edge. Default is 'reflect'. + imlib : str, optional + See the documentation of the ``vip_hci.preproc.frame_shift`` function. + interpolation : str, optional + See the documentation of the ``vip_hci.preproc.frame_shift`` function. debug : bool, optional If True debug information is printed and plotted (fit and residuals, intersections and shifts). This has to be used carefully as it can @@ -625,8 +634,9 @@ def cube_recenter_satspots(array, xy, subi_size=19, sigfactor=6, plot=True, res = frame_center_satspots(array[i], final_xy[i], debug=debug, shift=True, subi_size=subi_size, sigfactor=sigfactor, fit_type=fit_type, - filter_freq=filter_freq, - verbose=False, border_mode=border_mode) + filter_freq=filter_freq, imlib=imlib, + interpolation=interpolation, verbose=False, + border_mode=border_mode) array_rec.append(res[0]) shift_y[i] = res[1] shift_x[i] = res[2] @@ -1630,8 +1640,8 @@ def cube_recenter_2dfit(array, xy=None, fwhm=4, subi_size=5, model='gauss', plt.xlabel('Pixels') plt.figure(figsize=vip_figsize) - plt.plot(y, 'o-', label='shifts in y', alpha=0.5) plt.plot(x, 'o-', label='shifts in x', alpha=0.5) + plt.plot(y, 'o-', label='shifts in y', alpha=0.5) plt.legend(loc='best') plt.grid('on', alpha=0.2) plt.ylabel('Pixels') diff --git a/src/vip_hci/psfsub/pca_fullfr.py b/src/vip_hci/psfsub/pca_fullfr.py index dc195f21..e76634c5 100644 --- a/src/vip_hci/psfsub/pca_fullfr.py +++ b/src/vip_hci/psfsub/pca_fullfr.py @@ -120,8 +120,7 @@ class PCA_Params: smooth: float = None smooth_first_pass: float = None mask_rdi: np.ndarray = None - ref_strategy: str = 'RDI' # TBD: expand keyword to replace 'adimsdi' - # {'RDI', 'ARDI', 'RSDI', 'ARSDI', 'ASDI', 'S+ADI', 'S+ARDI'} + ref_strategy: str = 'RDI' # {'RDI', 'ARDI', 'RSDI', 'ARSDI'} check_memory: bool = True batch: Union[int, float] = None nproc: int = 1 @@ -130,6 +129,7 @@ class PCA_Params: weights: np.ndarray = None left_eigv: bool = False min_frames_pca: int = 10 + max_frames_pca: int = None cube_sig: np.ndarray = None med_of_npcs: bool = False @@ -373,9 +373,13 @@ def pca(*all_args: List, **all_kwargs: dict): [full_output=True, adimsdi='single'] Residuals cube (of the big cube with channels and time processed together). Valid for ADI+mSDI (4D) cubes (when ``scale_list`` is provided) + cube_desc_residuals : numpy ndarray + [full_output=True, adimsdi='single'] Residuals cube after de-scaling the + spectral frames to their original scale. Valid for ADI+mSDI (4D) (when + ``scale_list`` is provided). cube_adi_residuals : numpy ndarray - [full_output=True, adimsdi='single'] Residuals cube (of the big cube - with channels and time processed together) after de-scaling the wls. + [full_output=True, adimsdi='single'] Residuals cube after de-scaling the + spectral frames to their original scale and collapsing the channels. Valid for ADI+mSDI (4D) (when ``scale_list`` is provided). ifs_adi_frames : numpy ndarray [full_output=True, 4D input cube, ``scale_list=None``] This is the cube @@ -478,8 +482,21 @@ def pca(*all_args: List, **all_kwargs: dict): # ADI + mSDI. Shape of cube: (n_channels, n_adi_frames, y, x) # isinstance(cube, np.ndarray) and cube.ndim == 4: if algo_params.scale_list is not None: + add_params = {"start_time": start_time} + if algo_params.cube_ref is not None: + if algo_params.cube_ref.ndim != 4: + msg = "Ref cube has wrong format for 4d input cube" + raise TypeError(msg) + if 'A' in algo_params.ref_strategy: # e.g. 'ARSDI' + add_params["ref_strategy"] = 'ARSDI' # uniformize + if algo_params.adimsdi == Adimsdi.SINGLE: + cube_ref = np.concatenate((algo_params.cube, + algo_params.cube_ref), axis=1) + add_params["cube_ref"] = cube_ref + else: + add_params["ref_strategy"] = 'RSDI' + if algo_params.adimsdi == Adimsdi.DOUBLE: - add_params = {"start_time": start_time} func_params = setup_parameters( params_obj=algo_params, fkt=_adimsdi_doublepca, **add_params ) @@ -489,7 +506,6 @@ def pca(*all_args: List, **all_kwargs: dict): ) residuals_cube_channels, residuals_cube_channels_, frame = res_pca elif algo_params.adimsdi == Adimsdi.SINGLE: - add_params = {"start_time": start_time} func_params = setup_parameters( params_obj=algo_params, fkt=_adimsdi_singlepca, **add_params ) @@ -498,7 +514,8 @@ def pca(*all_args: List, **all_kwargs: dict): **rot_options, ) if np.isscalar(algo_params.ncomp): - cube_allfr_residuals, cube_adi_residuals, frame = res_pca + cube_allfr_residuals, cube_desc_residuals = res_pca[:2] + cube_adi_residuals, frame = res_pca[2:] elif isinstance(algo_params.ncomp, (tuple, list)): if algo_params.source_xy is None: if algo_params.full_output: @@ -700,7 +717,8 @@ def pca(*all_args: List, **all_kwargs: dict): # ADI+mSDI single-pass PCA if np.isscalar(algo_params.ncomp): if algo_params.full_output: - return frame, cube_allfr_residuals, cube_adi_residuals + return (frame, cube_allfr_residuals, cube_desc_residuals, + cube_adi_residuals) else: return frame # ADI+mSDI single-pass PCA grid @@ -791,6 +809,7 @@ def _adi_rdi_pca( cube_sig=None, left_eigv=False, min_frames_pca=10, + max_frames_pca=None, smooth=None, **rot_options, ): @@ -847,6 +866,9 @@ def _adi_rdi_pca( "Number of PCs too high (max PCs={}), using {} PCs " "instead.".format(nref, ncomp) ) + elif ncomp <= 0: + msg = "Number of PCs too low. It should be > 0." + raise ValueError(msg) if mask_rdi is None: if source_xy is None: residuals_result = _project_subtract( @@ -885,9 +907,16 @@ def _adi_rdi_pca( x1, y1 = source_xy ann_center = dist(yc, xc, y1, x1) pa_thr = _compute_pa_thresh(ann_center, fwhm, delta_rot) + max_fr = max_frames_pca + if max_frames_pca is not None: + truncate = True + else: + truncate = False for frame in range(n): - ind = _find_indices_adi(angle_list, frame, pa_thr) + ind = _find_indices_adi(angle_list, frame, pa_thr, + truncate=truncate, + max_frames=max_fr) res_result = _project_subtract( cube, @@ -995,6 +1024,7 @@ def _adi_rdi_pca( def _adimsdi_singlepca( cube, + cube_ref, angle_list, scale_list, ncomp, @@ -1018,6 +1048,7 @@ def _adimsdi_singlepca( weights=None, left_eigv=False, min_frames_pca=10, + ref_strategy='RSDI', **rot_options, ): """Handle the full-frame ADI+mSDI single PCA post-processing.""" @@ -1036,7 +1067,7 @@ def _adimsdi_singlepca( if not scale_list.shape[0] == z: raise ValueError("`scale_list` has wrong length") - scale_list = check_scal_vector(scale_list) + # scale_list = check_scal_vector(scale_list) big_cube = [] if verbose: @@ -1051,12 +1082,33 @@ def _adimsdi_singlepca( big_cube = np.array(big_cube) big_cube = big_cube.reshape(z * n, big_cube.shape[2], big_cube.shape[3]) + # Do the same with reference cube if provided + if cube_ref is not None: + z, nr, y_in, x_in = cube_ref.shape + big_cube_ref = [] + if verbose: + msg = "Rescaling the spectral channels of the reference cube.." + print(msg) + for i in Progressbar(range(nr), verbose=verbose): + cube_resc = scwave(cube_ref[:, i, :, :], scale_list, imlib=imlib2, + interpolation=interpolation)[0] + if crop_ifs: + cube_resc = cube_crop_frames(cube_resc, size=y_in, + verbose=False) + big_cube_ref.append(cube_resc) + + big_cube_ref = np.array(big_cube_ref) + big_cube_ref = big_cube_ref.reshape(z * nr, big_cube_ref.shape[2], + big_cube_ref.shape[3]) + else: + big_cube_ref = None + if verbose: timing(start_time) print("{} total frames".format(n * z)) print("Performing single-pass PCA") - if isinstance(ncomp, (int, float)): + if np.isscalar(ncomp): # When ncomp is a int and batch is not None, incremental ADI-PCA is run if batch is not None: res_cube = pca_incremental( @@ -1078,7 +1130,7 @@ def _adimsdi_singlepca( else: res_cube = _project_subtract( big_cube, - None, + big_cube_ref, ncomp, scaling, mask_center_px, @@ -1103,11 +1155,13 @@ def _adimsdi_singlepca( idx_ini = ifs_collapse_range[0] idx_fin = ifs_collapse_range[1] + cube_desc_residuals = np.zeros_like(cube) + for i in Progressbar(range(n), verbose=verbose): - frame_i = scwave( + res_i = scwave( res_cube[i * z + idx_ini:i * z + idx_fin], scale_list[idx_ini:idx_fin], - full_output=False, + full_output=True, inverse=True, y_in=y_in, x_in=x_in, @@ -1115,7 +1169,8 @@ def _adimsdi_singlepca( interpolation=interpolation, collapse=collapse_ifs, ) - resadi_cube[i] = frame_i + cube_desc_residuals[:, i] = res_i[0] + resadi_cube[i] = res_i[1] if verbose: print("De-rotating and combining residuals") @@ -1133,7 +1188,8 @@ def _adimsdi_singlepca( frame = cube_collapse(der_res, mode=collapse, w=weights) cube_allfr_residuals = res_cube cube_adi_residuals = resadi_cube - return cube_allfr_residuals, cube_adi_residuals, frame + return (cube_allfr_residuals, cube_desc_residuals, cube_adi_residuals, + frame) # When ncomp is a tuple, pca_grid is called elif isinstance(ncomp, (tuple, list)): @@ -1174,6 +1230,7 @@ def _adimsdi_singlepca( def _adimsdi_doublepca( cube, + cube_ref, angle_list, scale_list, ncomp, @@ -1195,14 +1252,22 @@ def _adimsdi_doublepca( delta_rot=None, fwhm=4, min_frames_pca=10, + max_frames_pca=None, mask_rdi=None, cube_sig=None, left_eigv=False, + ref_strategy='RSDI', **rot_options, ): """Handle the full-frame ADI+mSDI double PCA post-processing.""" z, n, y_in, x_in = cube.shape + if cube_ref is not None: + cube = np.concatenate((cube, cube_ref), axis=1) + nr = cube_ref.shape[1] + else: + nr = 0 + global ARRAY ARRAY = cube # to be passed to _adimsdi_doublepca_ifs @@ -1226,7 +1291,7 @@ def _adimsdi_doublepca( raise ValueError("Scaling factors vector is not 1d") if not scale_list.shape[0] == cube.shape[0]: raise ValueError("Scaling factors vector has wrong length") - scale_list = check_scal_vector(scale_list) + # scale_list = check_scal_vector(scale_list) if type(scaling) is not tuple: scaling = (scaling, scaling) @@ -1246,7 +1311,7 @@ def _adimsdi_doublepca( res = pool_map( nproc, _adimsdi_doublepca_ifs, - iterable(range(n)), + iterable(range(n+nr)), ncomp_ifs, scale_list, scaling[0], @@ -1277,7 +1342,7 @@ def _adimsdi_doublepca( print("{} ADI frames".format(n)) print("De-rotating and combining frames (skipping PCA)") residuals_cube_channels_ = cube_derotate( - res_cube_channels, + res_cube_channels[:n], angle_list, nproc=nproc, imlib=imlib, @@ -1289,27 +1354,43 @@ def _adimsdi_doublepca( if verbose: timing(start_time) else: - if ncomp_adi > n: - ncomp_adi = n + if ncomp_adi > n+nr: + ncomp_adi = n+nr msg = "Number of PCs too high, using maximum of {} PCs instead" print(msg.format(n)) if verbose: print("{} ADI frames".format(n)) + if nr: + print("+ {} reference frames".format(nr)) print("Second PCA stage exploiting rotational variability") if source_xy is None: - res_ifs_adi = _project_subtract( - res_cube_channels, - None, - ncomp_adi, - scaling[1], - mask_center_px, - svd_mode, - verbose, - False, - cube_sig=cube_sig, - left_eigv=left_eigv, - ) + if 'A' in ref_strategy or cube_ref is None: # e.g. 'ARSDI' + res_ifs_adi = _project_subtract( + res_cube_channels, + None, + ncomp_adi, + scaling[1], + mask_center_px, + svd_mode, + verbose, + False, + cube_sig=cube_sig, + left_eigv=left_eigv, + ) + else: # 'RSDI' + res_ifs_adi = _project_subtract( + res_cube_channels[:n], + res_cube_channels[n:], # ref cube + ncomp_adi, + scaling[1], + mask_center_px, + svd_mode, + verbose, + False, + cube_sig=cube_sig, + left_eigv=left_eigv, + ) if verbose: timing(start_time) # A rotation threshold is applied @@ -1324,12 +1405,19 @@ def _adimsdi_doublepca( pa_thr = _compute_pa_thresh(ann_center, fwhm, delta_rot) res_ifs_adi = np.zeros_like(res_cube_channels) + max_fr = max_frames_pca + if max_frames_pca is not None: + truncate = True + else: + truncate = False + for frame in range(n): - ind = _find_indices_adi(angle_list, frame, pa_thr) + ind = _find_indices_adi(angle_list, frame, pa_thr, + truncate=truncate, max_frames=max_fr) res_result = _project_subtract( - res_cube_channels, - None, + res_cube_channels[:n], + res_cube_channels[n:], ncomp_adi, scaling[1], mask_center_px, @@ -1346,15 +1434,10 @@ def _adimsdi_doublepca( if verbose: print("De-rotating and combining residuals") - der_res = cube_derotate( - res_ifs_adi, - angle_list, - nproc=nproc, - imlib=imlib, - interpolation=interpolation, - **rot_options, - ) - residuals_cube_channels_ = der_res + residuals_cube_channels_ = cube_derotate(res_ifs_adi[:n], angle_list, + nproc=nproc, imlib=imlib, + interpolation=interpolation, + **rot_options) frame = cube_collapse(residuals_cube_channels_, mode=collapse, w=weights) if verbose: diff --git a/src/vip_hci/psfsub/pca_local.py b/src/vip_hci/psfsub/pca_local.py index 39c3c49a..9293eec5 100644 --- a/src/vip_hci/psfsub/pca_local.py +++ b/src/vip_hci/psfsub/pca_local.py @@ -334,6 +334,10 @@ def pca_annular(*all_args: List, **all_kwargs: dict): global ARRAY ARRAY = algo_params.cube + if algo_params.cube_ref is not None: + global ARRAY_REF + ARRAY_REF = algo_params.cube_ref + z, n, y_in, x_in = algo_params.cube.shape algo_params.fwhm = int(np.round(np.mean(algo_params.fwhm))) n_annuli = int((y_in / 2 - algo_params.radius_int) / algo_params.asize) @@ -399,13 +403,44 @@ def pca_annular(*all_args: List, **all_kwargs: dict): ) else: + if algo_params.cube_ref is not None: + if algo_params.verbose: + print("First PCA subtraction (spectral) on REF cube") + print("{} spectral channels per IFS frame".format(z)) + print( + "N annuli = {}, mean FWHM = {:.3f}".format( + n_annuli, algo_params.fwhm) + ) + # apply the same first pass to cube ref + nr = algo_params.cube_ref.shape[0] + add_params['do_ref'] = True + add_params["fr"] = iterable(range(nr)) + func_params_ref = setup_parameters( + params_obj=algo_params, fkt=_pca_sdi_fr, as_list=True, + **add_params + ) + res = pool_map( + algo_params.nproc, + _pca_sdi_fr, + verbose=algo_params.verbose, + *func_params_ref, + ) + residuals_cube_channels_ref = np.array(res) + + # Exploiting rotational variability + if algo_params.verbose: + timing(start_time) + print("{} REF frames".format(n)) + else: + residuals_cube_channels_ref = None + if algo_params.verbose: print("Second PCA subtraction exploiting angular variability") add_params = { "cube": residuals_cube_channels, "ncomp": ncomp2, - "cube_ref": None, + "cube_ref": residuals_cube_channels_ref, } func_params = setup_parameters( @@ -449,15 +484,26 @@ def _pca_sdi_fr( collapse, ifs_collapse_range, theta_init, + do_ref=False ): """Optimized PCA subtraction on a multi-spectral frame (IFS data).""" - z, n, y_in, x_in = ARRAY.shape - scale_list = check_scal_vector(scal) - # rescaled cube, aligning speckles - multispec_fr = scwave( - ARRAY[:, fr, :, :], scale_list, imlib=imlib, interpolation=interpolation - )[0] + + if do_ref: # do it on REF cube + z, n, y_in, x_in = ARRAY_REF.shape + # rescaled cube, aligning speckles + multispec_fr = scwave( + ARRAY_REF[:, fr, :, :], scale_list, imlib=imlib, + interpolation=interpolation + )[0] + else: # do it on SCI cube + z, n, y_in, x_in = ARRAY.shape + + # rescaled cube, aligning speckles + multispec_fr = scwave( + ARRAY[:, fr, :, :], scale_list, imlib=imlib, + interpolation=interpolation + )[0] # Exploiting spectral variability (radial movement) fwhm = int(np.round(np.mean(fwhm))) diff --git a/src/vip_hci/stats/clip_sigma.py b/src/vip_hci/stats/clip_sigma.py index 59f2dd9a..d5afaf77 100644 --- a/src/vip_hci/stats/clip_sigma.py +++ b/src/vip_hci/stats/clip_sigma.py @@ -191,8 +191,8 @@ def _sigma_filter(frame_tmp, bpix_map, neighbor_box=3, min_neighbors=3, min_neighbors=3, verbose=False) -def clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori=None, - out_good=False, neighbor=False, num_neighbor=3, mad=False, +def clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori=None, + out_good=False, neighbor=False, num_neighbor=3, mad=False, half_res_y=False): """Sigma clipping for detecting outlying values in 2d array. If the parameter 'neighbor' is True the clipping can be performed in a local patch @@ -211,7 +211,7 @@ def clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori=None, on top of those ones. They are not considered as good neighbours during sigma-clipping. out_good : bool, optional - For choosing different outputs. True to return indices of good pixels, + For choosing different outputs. True to return indices of good pixels, False for indices of bad pixels. neighbor : bool, optional For clipping over the median of the contiguous pixels. @@ -239,11 +239,11 @@ def clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori=None, """ - def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, + def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, out_good=False, neighbor=False, num_neighbor=3, mad=False, half_res_y=False): """Sigma clipping for detecting outlying values in 2d array. If the - parameter 'neighbor' is True the clipping can be performed in a local + parameter 'neighbor' is True the clipping can be performed in a local patch around each pixel, whose size depends on 'neighbor' parameter. Parameters @@ -255,8 +255,8 @@ def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, upper_sigma : float Value for sigma, upper boundary. bpm_mask_ori : 2d numpy ndarray or None - Known (e.g. static) bad pixel map. Additional bad pixels will be - identified, on top of these ones. They are not considered as good + Known (e.g. static) bad pixel map. Additional bad pixels will be + identified, on top of these ones. They are not considered as good neighbours during sigma-clipping. out_good : bool, optional For choosing different outputs. @@ -288,7 +288,7 @@ def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, gpm_ori = np.ones(array.shape) else: gpm_ori = np.ones(array.shape) - bpm_mask_ori - + ny, nx = array.shape bpm = np.ones(array.shape) gpm = np.zeros(array.shape) @@ -373,15 +373,15 @@ def _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, else: bad = np.where(bpm) return bad - + if no_numba: - return _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, - out_good=out_good, neighbor=neighbor, - num_neighbor=num_neighbor, mad=mad, + return _clip_array(array, lower_sigma, upper_sigma, bpm_mask_ori, + out_good=out_good, neighbor=neighbor, + num_neighbor=num_neighbor, mad=mad, half_res_y=half_res_y) else: _clip_array_numba = njit(_clip_array) return _clip_array_numba(array, lower_sigma, upper_sigma, bpm_mask_ori, out_good=out_good, neighbor=neighbor, num_neighbor=num_neighbor, mad=mad, - half_res_y=half_res_y) \ No newline at end of file + half_res_y=half_res_y) diff --git a/src/vip_hci/var/iuwt.py b/src/vip_hci/var/iuwt.py index b7db0f6a..a19bcc17 100644 --- a/src/vip_hci/var/iuwt.py +++ b/src/vip_hci/var/iuwt.py @@ -14,10 +14,10 @@ def iuwt_decomposition(in1, scale_count, scale_adjust=0, mode='ser', core_count=2, store_smoothed=False): """ - This function serves as a handler for the different implementations of the - IUWT decomposition. It allows the different methods to be used almost + This function serves as a handler for the different implementations of the + IUWT decomposition. It allows the different methods to be used almost interchangeably. - + The code was taken from [KEN15]_ and is detailed in [DAB15]_. INPUTS: diff --git a/tests/pre_3_10/test_fm_negfc_3d.py b/tests/pre_3_10/test_fm_negfc_3d.py index 1d457efe..9006a8b0 100644 --- a/tests/pre_3_10/test_fm_negfc_3d.py +++ b/tests/pre_3_10/test_fm_negfc_3d.py @@ -122,7 +122,7 @@ def test_algos( algo_options["radius_int"] = res0[0][0] - 2 * ds.fwhm algo_options["asize"] = 4 * ds.fwhm algo_options["delta_rot"] = 1 - if pca_algo == pca: + if pca_algo == median_sub: # just test it once because very slow sp_unc = speckle_noise_uncertainty( cube_emp,