diff --git a/docs/source/conf.py b/docs/source/conf.py index 5d265681..11be4169 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -53,6 +53,8 @@ "dollarmath" ] +suppress_warnings = ["mystnb.unknown_mime_type"] + # Add any paths that contain templates here, relative to this directory. templates_path = []#['./_templates'] diff --git a/docs/source/tutorials/01A_quickstart.ipynb b/docs/source/tutorials/01A_quickstart.ipynb index 849a0219..fa47b1cc 100644 --- a/docs/source/tutorials/01A_quickstart.ipynb +++ b/docs/source/tutorials/01A_quickstart.ipynb @@ -105,7 +105,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "```{warning} \n", + "```{important} \n", "Your version of VIP should be >=1.0.3. If not, close the notebook and run for example `pip install vip_hci --upgrade` in your terminal.\n", "```" ] @@ -198,7 +198,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To inspect or retrieve the data, the easiest way is to provide their local path to the relevant functions in VIP - as illustrated below.\n", + "To inspect or retrieve the data, the easiest way is to provide their local path to the relevant functions in VIP - as illustrated with the commented code lines below.\n", "\n", "Here we use the `download_file` utility of `astropy` with the `cache` option - which equivalently stores the local path to the file after downloading the data from an online repository:" ] @@ -260,7 +260,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "```{warning} \n", + "```{important} \n", "The *convention* in VIP is to always consider you are providing **derotation angles** as input, not **parallactic angles**. Derotation angles are the angles to be applied to the observed images (with positive offsets measured counter-clockwise) in order to align them with north up. These are essentially the *opposite* values of the parallactic angles plus some instrument-specific offsets. The additional offsets can correspond to any internal reflections affecting the orientation of the field of view, or to a correction for the true north orientation which may drift over time and requires dedicated instrumental calibrations. These additional offsets should be checked in the manual of the instrument you are using.\n", "```" ] @@ -1351,7 +1351,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "```{warning} \n", + "```{important} \n", "VIP's convention regarding star centering is:\n", "\n", "- for **odd** number of pixels along the x and y dimensions: the star centroid should correspond to the central pixel of the image;\n", diff --git a/docs/source/tutorials/01B_quickstart_with_objects.ipynb b/docs/source/tutorials/01B_quickstart_with_objects.ipynb index fc6dd4d3..e9d5ddca 100644 --- a/docs/source/tutorials/01B_quickstart_with_objects.ipynb +++ b/docs/source/tutorials/01B_quickstart_with_objects.ipynb @@ -22,11 +22,11 @@ "source": [ "**Table of contents**\n", "\n", - "* [1.1. Loading ADI data](#1.1.-Loading-ADI-data)\n", - "* [1.2. Model PSF subtraction for ADI](#1.2.-Model-PSF-subtraction-for-ADI)\n", - " - [1.2.1. median-ADI](#1.2.1.-median-ADI)\n", - " - [1.2.2. Full-frame PCA](#1.2.2.-Full-frame-PCA) \n", - "* [1.3. Is the point source significant?](#1.3.-Is-it-significant?)" + "* [1.1. Loading ADI data](1B.1.-Loading-ADI-data)\n", + "* [1.2. Model PSF subtraction for ADI](1B.2.-Model-PSF-subtraction-for-ADI)\n", + " - [1.2.1. median-ADI](1B.2.1.-median-ADI)\n", + " - [1.2.2. Full-frame PCA](1B.2.2.-Full-frame-PCA) \n", + "* [1.3. Is the point source significant?](1B.3.-Is-it-significant?)" ] }, { @@ -108,6 +108,15 @@ " raise ValueError(msg)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{important}\n", + "Your version of VIP should be >=1.4.2 to run this notebook. If not, close the notebook and run for example `pip install vip_hci --upgrade` in your terminal. Also make sure the Python version used is 3.10 or higher. \n", + "```" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -119,6 +128,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "(1B.1.-Loading-ADI-data)=\n", "## 1.1. Loading ADI data" ] }, @@ -126,8 +136,32 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In order to demonstrate the capabilities of ``VIP``, you can find in the 'dataset' folder of the `VIP_extras` repository: a toy Angular Differential Imaging (ADI) datacube of images and a point spread function (PSF) associated to the NACO instrument of the Very Large Telescope (VLT). This is an L'-band VLT/NACO dataset of beta Pictoris published in [Absil et al. (2013)](https://ui.adsabs.harvard.edu/abs/2013A%26A...559L..12A/abstract), obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells.\n", + "In order to demonstrate the capabilities of ``VIP``, you can find in the 'dataset' folder of the `VIP_extras` repository: a toy Angular Differential Imaging (ADI) datacube of images and a point spread function (PSF) associated to the NACO instrument of the Very Large Telescope (VLT). This is an L'-band VLT/NACO dataset of beta Pictoris published in [Absil et al. (2013)](https://ui.adsabs.harvard.edu/abs/2013A%26A...559L..12A/abstract), obtained using the Annular Groove Phase Mask (AGPM) Vortex coronagraph. The sequence has been heavily sub-sampled temporarily to make it smaller. The frames were also cropped to the central 101x101 area. In case you want to plug-in your cube just change the path of the following cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.utils.data import download_file\n", + "\n", + "url_d = \"https://github.com/vortex-exoplanet/VIP_extras/raw/master/datasets\"\n", + "f1 = download_file(\"{}/naco_betapic_cube_cen.fits\".format(url_d), cache=True)\n", + "f2 = download_file(\"{}/naco_betapic_psf.fits\".format(url_d), cache=True)\n", + "f3 = download_file(\"{}/naco_betapic_derot_angles.fits\".format(url_d), cache=True)\n", "\n", + "# alternatively, for local files simply provide their full or relative path. E.g.:\n", + "#f1 = '../datasets/naco_betapic_cube_cen.fits'\n", + "#f2 = '../datasets/naco_betapic_psf.fits'\n", + "#f3 = '../datasets/naco_betapic_derot_angles.fits'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "Let's first inspect the fits file:" ] }, @@ -151,7 +185,7 @@ "source": [ "from vip_hci.fits import info_fits\n", "\n", - "info_fits('../datasets/naco_betapic_cube_cen.fits')" + "info_fits(f1)" ] }, { @@ -160,20 +194,23 @@ "source": [ "The fits file contains an ADI datacube of 61 images 101x101 px across.\n", "\n", - "The parallactic angle vector for this observation is saved in the `naco_betapic_pa.fits` file. \n", - "\n", - "Now we load all the data in memory with the `Dataset` object:" + "The parallactic angle vector for this observation is saved in the `naco_betapic_pa.fits` file. " ] }, { - "cell_type": "code", - "execution_count": 4, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "psfname = '../datasets/naco_betapic_psf.fits'\n", - "cubename = '../datasets/naco_betapic_cube_cen.fits'\n", - "angname = '../datasets/naco_betapic_pa.fits'" + "```{important} \n", + "The *convention* in VIP is to always consider you are providing **derotation angles** as input, not **parallactic angles**. Derotation angles are the angles to be applied to the observed images (with positive offsets measured counter-clockwise) in order to align them with north up. These are essentially the *opposite* values of the parallactic angles plus some instrument-specific offsets. The additional offsets can correspond to any internal reflections affecting the orientation of the field of view, or to a correction for the true north orientation which may drift over time and requires dedicated instrumental calibrations. These additional offsets should be checked in the manual of the instrument you are using.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we load all the data in memory with the `Dataset` object:" ] }, { @@ -193,7 +230,7 @@ ], "source": [ "from vip_hci.objects import Dataset\n", - "betapic = Dataset(cube=cubename, psf=psfname, angles=angname)" + "betapic = Dataset(cube=f1, psf=f2, angles=f3)" ] }, { @@ -1156,12 +1193,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "---\n", - "**NOTE**\n", - "\n", + "```{note}\n", "The `normalize_psf` function performs internally a fine centering of the PSF. The input PSF should nevertheless already be centered within a few pixel accuracy for the fine centering to work.\n", - "\n", - "---" + "```" ] }, { @@ -1223,7 +1257,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 1.2. Model PSF subtraction for ADI" + "(1B.2.-PSF-subtraction)=\n", + "## 1.2. PSF subtraction" ] }, { @@ -1237,15 +1272,19 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "**IMPORTANT** - VIP's convention regarding centering is:\n", - "- for **odd** number of pixels along the x and y directions: the star is placed on the central pixel;\n", - "- for **even** number of pixels: the star is placed on coordinates (dim/2,dim/2) in 0-based indexing." + "```{important} \n", + "VIP's convention regarding star centering is:\n", + "\n", + "- for **odd** number of pixels along the x and y dimensions: the star centroid should correspond to the central pixel of the image;\n", + "- for **even** number of pixels: the star centroid should be located on coordinates (dim/2,dim/2) in 0-based indexing (i.e. on the top-right pixel of the 4 central pixels).\n", + "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "(1B.2.1.-median-ADI)=\n", "### 1.2.1. median-ADI" ] }, @@ -1365,6 +1404,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "```{note}\n", + "As per usual convention, north is up and east is left in all images shown in the tutorial. Note that for the images of some instruments, te x-axis may be flipped (i.e. east would be right instead of left). This is for example the case for Gemini/NICI or VLT/ERIS-NIX. Double-check your instrument manual.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(1B.2.2.-Full-frame-PCA)=\n", "### 1.2.2. Full-frame PCA" ] }, @@ -1470,6 +1519,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "(1B.3.-Is-it-significant?)=\n", "## 1.3. Is it significant?" ] }, @@ -1575,6 +1625,13 @@ "toc_position": {}, "toc_section_display": true, "toc_window_display": true + }, + "widgets": { + "application/vnd.jupyter.widget-state+json": { + "state": {}, + "version_major": 2, + "version_minor": 0 + } } }, "nbformat": 4, diff --git a/vip_hci/fm/fakecomp.py b/vip_hci/fm/fakecomp.py index b90e0666..9eb28492 100644 --- a/vip_hci/fm/fakecomp.py +++ b/vip_hci/fm/fakecomp.py @@ -793,7 +793,7 @@ def psf_norm_2d(psf, fwhm, threshold, mask_core, full_output, verbose): restemp = psf_norm_2d(array[fr], fwhm[fr], threshold, mask_core, True, False) array_out.append(restemp[0]) - fwhm_flux[fr] = restemp[1] + fwhm_flux[fr] = restemp[1][0] array_out = np.array(array_out) if verbose: diff --git a/vip_hci/fm/negfc_fmerit.py b/vip_hci/fm/negfc_fmerit.py index 08b0c6e2..88278007 100644 --- a/vip_hci/fm/negfc_fmerit.py +++ b/vip_hci/fm/negfc_fmerit.py @@ -305,7 +305,7 @@ def chisquare( msg += "calculation. Consider larger input images." raise ValueError(msg) - subim = frame_crop(frpca, crop_sz, cenxy=(xround, yround), + subim = frame_crop(frpca, crop_sz, xy=(xround, yround), force=True, verbose=False) H = hessian(subim) dets = np.zeros([ndet, ndet]) diff --git a/vip_hci/invprob/fmmf.py b/vip_hci/invprob/fmmf.py index 777a40f9..bc17eb52 100644 --- a/vip_hci/invprob/fmmf.py +++ b/vip_hci/invprob/fmmf.py @@ -460,7 +460,7 @@ def _snr_contrast_esti( psfm = frame_crop( psfm_temp[j], crop, - cenxy=[int(psfm_temp.shape[-1] / 2), + xy=[int(psfm_temp.shape[-1] / 2), int(psfm_temp.shape[-1] / 2)], verbose=False, ) @@ -468,7 +468,7 @@ def _snr_contrast_esti( num.append( np.multiply( frame_crop( - mcube[j], crop, cenxy=[ + mcube[j], crop, xy=[ poscentx, poscenty], verbose=False ), psfm, diff --git a/vip_hci/metrics/contrcurve.py b/vip_hci/metrics/contrcurve.py index 8aa5a350..953cbb6b 100644 --- a/vip_hci/metrics/contrcurve.py +++ b/vip_hci/metrics/contrcurve.py @@ -703,12 +703,12 @@ def throughput( raise TypeError(msg) if psf_template.ndim != 3: raise TypeError("Template PSF is not a frame, 3d array") - if "scale_list" not in algo_dict: - raise ValueError("Vector of wavelength not found") - else: + if "scale_list" in algo_dict: + # raise ValueError("Vector of wavelength not found") + # else: if algo_dict["scale_list"].shape[0] != array.shape[0]: raise TypeError("Input wavelength vector has wrong length") - if isinstance(fwhm, float) or isinstance(fwhm, int): + if np.isscalar(fwhm): maxfcsep = int((array.shape[2] / 2.0) / fwhm) - 1 else: maxfcsep = int((array.shape[2] / 2.0) / np.amin(fwhm)) - 1 diff --git a/vip_hci/metrics/snr_source.py b/vip_hci/metrics/snr_source.py index 747824b7..34c61cd1 100644 --- a/vip_hci/metrics/snr_source.py +++ b/vip_hci/metrics/snr_source.py @@ -452,7 +452,8 @@ def snr(array, source_xy, fwhm, full_output=False, array2=None, use2alone=False, return snr_vale -def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True): +def significance(snr, rad, fwhm, n_ap=None, student_to_gauss=True, + verbose=True): """Convert a S/N ratio (measured as in [MAW14]_) into a Gaussian\ significance (n-sigma) with equivalent false alarm probability for a\ point-source detection measured at a given separation, or the opposite. @@ -467,6 +468,11 @@ def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True): to each snr measurement. fwhm : float Full Width Half Maximum of the PSF. + n_ap : int, optional + Number of independent samples available. If provided, will take + precedence over using input rad and FWHM. This can be useful when the + limited azimuthal coverage is available (e.g. when using APP or 4QPM + coronagraphs). student_to_gauss : bool, optional Whether the conversion is from Student SNR to Gaussian significance (True), or the opposite (False). @@ -477,9 +483,12 @@ def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True): Equivalent Gaussian significance [student_to_gauss=True] or equivalent Student S/N ratio [student_to_gauss=False]. """ + if n_ap is None: + n_ap = (rad/fwhm)*2*np.pi-2 + if student_to_gauss: - sig = norm.ppf(t.cdf(snr, (rad/fwhm)*2*np.pi-2)) - if t.cdf(snr, (rad/fwhm)*2*np.pi-2) == 1.0: + sig = norm.ppf(t.cdf(snr, n_ap)) + if t.cdf(snr, n_ap) == 1.0: print("Warning high S/N! cdf>0.9999999999999999 is rounded to 1") msg = "Returning 8.2 sigma, but quote significance > 8.2 sigma." print(msg) @@ -490,7 +499,7 @@ def significance(snr, rad, fwhm, student_to_gauss=True, verbose=True): msg += r"Gaussian false alarm probability." print(msg.format(rad, rad/fwhm, snr, sig)) else: - sig = t.ppf(norm.cdf(snr), (rad/fwhm)*2*np.pi-2) + sig = t.ppf(norm.cdf(snr), n_ap) if verbose: msg = r"At a separation of {:.1f} px ({:.1f} FWHM), a {:.1f}-sigma " msg += r"detection in terms of Gaussian false alarm probability " diff --git a/vip_hci/preproc/cosmetics.py b/vip_hci/preproc/cosmetics.py index ae819690..65dcf81f 100644 --- a/vip_hci/preproc/cosmetics.py +++ b/vip_hci/preproc/cosmetics.py @@ -4,7 +4,6 @@ Also functions for cropping cubes. """ - __author__ = "Carlos Alberto Gomez Gonzalez, V. Christiaens" __all__ = [ "cube_crop_frames", @@ -15,16 +14,12 @@ "approx_stellar_position", ] - from multiprocessing import cpu_count import numpy as np from astropy.stats import sigma_clipped_stats -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 set_start_method shared_mem = True try: @@ -79,9 +74,8 @@ def cube_crop_frames( cenx, ceny = xy else: ceny, cenx = frame_center(temp_fr) - _, y0, x0 = get_square( - temp_fr, size, y=ceny, x=cenx, position=True, force=force, verbose=verbose - ) + _, y0, x0 = get_square(temp_fr, size, y=ceny, x=cenx, position=True, + force=force, verbose=verbose) if not force: if temp_fr.shape[0] % 2 == 0: @@ -107,7 +101,7 @@ def cube_crop_frames( return array_out -def frame_crop(array, size, cenxy=None, force=False, verbose=True): +def frame_crop(array, size, xy=None, force=False, verbose=True): """Crops a square frame (2d array). Uses the ``get_square`` function. Parameters @@ -116,7 +110,7 @@ def frame_crop(array, size, cenxy=None, force=False, verbose=True): Input frame. size : int, odd Size of the subframe. - cenxy : tuple, optional + xy : tuple, optional Coordinates of the center of the subframe. force : bool, optional Size and the size of the 2d array must be both even or odd. With @@ -136,11 +130,12 @@ def frame_crop(array, size, cenxy=None, force=False, verbose=True): if array.ndim != 2: raise TypeError("`Array` is not a frame or 2d array") - if not cenxy: + if not xy: ceny, cenx = frame_center(array) else: - cenx, ceny = cenxy - array_view = get_square(array, size, ceny, cenx, force=force, verbose=verbose) + cenx, ceny = xy + array_view = get_square(array, size, ceny, cenx, force=force, + verbose=verbose) if verbose: print("New shape: {}".format(array_view.shape)) @@ -267,16 +262,16 @@ def cube_drop_frames(array, n, m, parallactic=None, verbose=True): raise TypeError("End index must be smaller than the # of frames") if array.ndim == 4: - array_view = array[:, n - 1 : m, :, :].copy() + array_view = array[:, n - 1: m, :, :].copy() elif array.ndim == 3: - array_view = array[n - 1 : m, :, :].copy() + array_view = array[n - 1: m, :, :].copy() else: raise ValueError("only 3D and 4D cubes are supported!") if parallactic is not None: if not parallactic.ndim == 1: raise ValueError("Parallactic angles vector has wrong shape") - parallactic = parallactic[n - 1 : m] + parallactic = parallactic[n - 1: m] if verbose: print("Cube successfully sliced") @@ -292,9 +287,8 @@ def cube_drop_frames(array, n, m, parallactic=None, verbose=True): def frame_remove_stripes(array): - """Removes unwanted stripe artifact in frames with non-perfect bias or sky - subtraction. Encountered this case on an LBT data cube. - """ + """Remove unwanted stripe artifact in frames with non-perfect bias or sky\ + subtraction. Encountered this case on an LBT data cube.""" lines = array[:50] lines = np.vstack((lines, array[-50:])) mean = lines.mean(axis=0) @@ -303,11 +297,10 @@ def frame_remove_stripes(array): return array -def cube_correct_nan( - cube, neighbor_box=3, min_neighbors=3, verbose=False, half_res_y=False, nproc=1 -): - """Sigma filtering of nan pixels in a whole frame or cube. Tested on - SINFONI data. +def cube_correct_nan(cube, neighbor_box=3, min_neighbors=3, verbose=False, + half_res_y=False, nproc=1): + """ + Sigma filter nan pixels in a whole frame or cube. Parameters ---------- @@ -334,8 +327,8 @@ def cube_correct_nan( ------- obj_tmp : numpy ndarray Output cube with corrected nan pixels in each frame - """ + """ obj_tmp = cube.copy() ndims = obj_tmp.ndim @@ -364,14 +357,15 @@ def cube_correct_nan( if nproc == 1 or not shared_mem: for zz in range(n_z): obj_tmp[zz], nnanpix = nan_corr_2d( - obj_tmp[zz], neighbor_box, min_neighbors, half_res_y, verbose, True + obj_tmp[zz], neighbor_box, min_neighbors, half_res_y, + verbose, True ) if verbose: 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 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, @@ -381,16 +375,17 @@ def cube_correct_nan( full_output=False, ) if verbose: - msg = "Correcting NaNs in multiprocessing using ADACS' approach..." + msg = "Correcting NaNs in multiprocessing using ADACS' approach" print(msg) # 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. + # creating array object similar to obj_tmp and occupying shared mem 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. + # 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 ): @@ -417,7 +412,8 @@ def _nan_corr_2d_mp(args): args = [] for j in range(n_z): args.append( - [j, obj_tmp[j], neighbor_box, min_neighbors, half_res_y, verbose] + [j, obj_tmp[j], neighbor_box, min_neighbors, half_res_y, + verbose] ) try: pool.map_async(_nan_corr_2d_mp, args, chunksize=1).get( @@ -429,9 +425,11 @@ def _nan_corr_2d_mp(args): 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) + # 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: @@ -487,9 +485,9 @@ def nan_corr_2d( def approx_stellar_position(cube, fwhm, return_test=False, verbose=False): - """Finds the approximate coordinates of the star, assuming it is the - brightest signal in the images. The algorithm can handle images dominated - by noise, since outliers are corrected based on the position of ths star in + """Find the approximate coordinates of the star, assuming it is the\ + brightest signal in the images. The algorithm can handle images dominated\ + by noise, since outliers are corrected based on the position of ths star in\ other channels. Parameters diff --git a/vip_hci/psfsub/pca_fullfr.py b/vip_hci/psfsub/pca_fullfr.py index 91695c60..62eacb7b 100644 --- a/vip_hci/psfsub/pca_fullfr.py +++ b/vip_hci/psfsub/pca_fullfr.py @@ -108,8 +108,6 @@ class PCA_Params: source_xy: Tuple[int] = None delta_rot: int = None fwhm: float = 4 - # strategy: str = 'ADI' # TBD: add a strategy keyword: 'ADI', 'RDI', 'ARDI', - # 'ASDI', 'SDI', 'S+ADI', 'ARSDI', 'RSDI' => replace 'adimsdi' adimsdi: Enum = Adimsdi.SINGLE crop_ifs: bool = True imlib: Enum = Imlib.VIPFFT @@ -119,6 +117,8 @@ class PCA_Params: collapse_ifs: Enum = Collapse.MEAN ifs_collapse_range: Union[str, Tuple[int]] = "all" mask_rdi: np.ndarray = None + ref_strategy: str = 'RDI' # TBD: expand keyword to replace 'adimsdi' + # {'RDI', 'ARDI', 'RSDI', 'ARSDI', 'ASDI', 'S+ADI', 'S+ARDI'} check_memory: bool = True batch: Union[int, float] = None nproc: int = 1 @@ -287,6 +287,11 @@ def pca(*all_args: List, **all_kwargs: dict): and boat regions, respectively, following the denominations in [REN23]_. If only one mask is provided, it will be used as the anchor, and the boat images will not be masked (i.e., full frames used). + ref_strategy: str, opt {'RDI', 'ARDI'} + [cube_ref is not None] Indicates the strategy to be adopted when a + reference cube is provided. By default, RDI is done - i.e. the science + images are not used in the PCA library. If set to 'ARDI', the PCA + library is made of both the science and reference images. check_memory : bool, optional If True, it checks that the input cube is smaller than the available system memory. @@ -473,7 +478,7 @@ def pca(*all_args: List, **all_kwargs: dict): **func_params, **rot_options, ) - if isinstance(algo_params.ncomp, (int, float)): + if np.isscalar(algo_params.ncomp): cube_allfr_residuals, cube_adi_residuals, frame = res_pca elif isinstance(algo_params.ncomp, (tuple, list)): if algo_params.source_xy is None: @@ -527,7 +532,17 @@ def pca(*all_args: List, **all_kwargs: dict): if algo_params.cube_ref[ch].ndim != 3: msg = "Ref cube has wrong format for 4d input cube" raise TypeError(msg) - add_params["cube_ref"] = algo_params.cube_ref[ch] + if algo_params.ref_strategy == 'RDI': + add_params["cube_ref"] = algo_params.cube_ref[ch] + elif algo_params.ref_strategy == 'ARDI': + cube_ref = np.concatenate((algo_params.cube[ch], + algo_params.cube_ref[ch])) + add_params["cube_ref"] = cube_ref + else: + msg = "ref_strategy argument not recognized." + msg += "Should be 'RDI' or 'ARDI'" + raise TypeError(msg) + func_params = setup_parameters( params_obj=algo_params, fkt=_adi_rdi_pca, **add_params @@ -601,11 +616,19 @@ def pca(*all_args: List, **all_kwargs: dict): "full_output": True, } - func_params = setup_parameters(params_obj=algo_params, - fkt=_adi_rdi_pca, **add_params) - if algo_params.cube_ref is not None and algo_params.batch is not None: raise ValueError("RDI not compatible with batch mode") + elif algo_params.cube_ref is not None: + if algo_params.ref_strategy == 'ARDI': + algo_params.cube_ref = np.concatenate((algo_params.cube, + algo_params.cube_ref)) + elif algo_params.ref_strategy != 'RDI': + msg = "ref_strategy argument not recognized." + msg += "Should be 'RDI' or 'ARDI'" + raise TypeError(msg) + + func_params = setup_parameters(params_obj=algo_params, + fkt=_adi_rdi_pca, **add_params) res_pca = _adi_rdi_pca(**func_params, **rot_options) @@ -652,7 +675,7 @@ def pca(*all_args: List, **all_kwargs: dict): elif algo_params.adimsdi == Adimsdi.SINGLE: # ADI+mSDI single-pass PCA - if isinstance(algo_params.ncomp, (float, int)): + if np.isscalar(algo_params.ncomp): if algo_params.full_output: return frame, cube_allfr_residuals, cube_adi_residuals else: @@ -781,7 +804,7 @@ def _adi_rdi_pca( "equal the number of frames in the cube" ) - if not isinstance(ncomp, (int, float, tuple, list)): + if not np.isscalar(ncomp) and not isinstance(ncomp, (tuple, list)): msg = "`ncomp` must be an int, float, tuple or list in the ADI case" raise TypeError(msg) diff --git a/vip_hci/psfsub/pca_local.py b/vip_hci/psfsub/pca_local.py index 4b19bcb8..b71c6fd2 100644 --- a/vip_hci/psfsub/pca_local.py +++ b/vip_hci/psfsub/pca_local.py @@ -38,11 +38,7 @@ @dataclass class PCA_ANNULAR_Params: - """ - Set of parameters for the annular PCA module. - - - """ + """Set of parameters for the annular PCA module.""" cube: np.ndarray = None angle_list: np.ndarray = None