Merge branch 'master' of https://github.com/vortex-exoplanet/VIP

README.rst

@@ -70,8 +70,10 @@ Tutorials, in the form of Jupyter notebooks, showcasing ``VIP``'s usage and
 other resources such as test datasets are available in the 
 ``VIP-extras`` `repository <https://github.com/vortex-exoplanet/VIP_extras>`_. 
 Alternatively, you can execute this repository on 
-`Binder <https://mybinder.org/v2/gh/vortex-exoplanet/VIP_extras/master>`_ (in the tutorials directory). The first (quick-start) notebook can be visualized online with
-`nbviewer <http://nbviewer.jupyter.org/github/vortex-exoplanet/VIP_extras/blob/master/tutorials/01_quickstart.ipynb>`_. 
+`Binder <https://mybinder.org/v2/gh/vortex-exoplanet/VIP_extras/master>`_ (in 
+the tutorials directory). The first (quick-start) notebook can be visualized 
+online with `nbviewer 
+<http://nbviewer.jupyter.org/github/vortex-exoplanet/VIP_extras/blob/master/tutorials/01_quickstart.ipynb>`_. 
 If you are new to the Jupyter notebook application check out the `beginner's guide
 <https://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/what_is_jupyter.html>`_.
 

requirements-dev.txt

@@ -7,4 +7,5 @@ codecov ~=2.0.15
 flake8 ~= 3.5.0
 flake8-bandit ~= 1.0.2
 flake8-docstrings ~= 1.3.0
-autopep8
+autopep8
+sphinx-rtd-theme

vip_hci/__init__.py

@@ -1,5 +1,6 @@
 __version__ = "1.3.2"
 
+
 from . import preproc
 from . import config
 from . import fits

vip_hci/fits/fits.py

@@ -72,7 +72,7 @@ def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False,
         data = np.array(data, dtype=precision)
 
         if header:
-            header = hdulist[0].header
+            header = hdulist[n].header
             if verbose:
                 print("Fits HDU-{} data and header successfully loaded. "
                       "Data shape: {}".format(n, data.shape))

vip_hci/fm/negfc_fmerit.py

@@ -11,7 +11,8 @@
 from hciplot import plot_frames
 from skimage.draw import disk
 from ..fm import cube_inject_companions
-from ..var import (frame_center, get_annular_wedge, cube_filter_highpass)
+from ..var import (frame_center, get_annular_wedge, cube_filter_highpass,
+                   get_annulus_segments)
 from ..psfsub import pca_annulus, pca_annular, pca
 from ..preproc import cube_crop_frames
 
@@ -357,7 +358,25 @@ def get_values_optimize(cube, angs, ncomp, annulus_width, aperture_radius,
 
     indices = disk((posy, posx), radius=aperture_radius*fwhm)
     yy, xx = indices
-
+
+    # also consider indices of the annulus for pca_annulus
+    if algo == pca_annulus:
+        fr_size = res.shape[-1]
+        inner_rad = r_guess-annulus_width/2
+        yy_a, xx_a = get_annulus_segments((fr_size, fr_size), inner_rad,
+                                          annulus_width, nsegm=1)[0]
+        ## only consider overlapping indices
+        yy_f = []
+        xx_f = []
+        for i in range(len(yy)):
+            ind_y = np.where(yy_a==yy[i])
+            for j in ind_y[0]:
+                if xx[i]==xx_a[j]:
+                    yy_f.append(yy[i])
+                    xx_f.append(xx[i])
+        yy = np.array(yy_f, dtype=int)
+        xx = np.array(xx_f, dtype=int)
+
     if collapse is None:
         values = res[:, yy, xx].ravel()
     else:

vip_hci/fm/negfc_speckle_noise.py

@@ -208,7 +208,7 @@ def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo,
         if len(mu_sigma) != 2:
             raise TypeError("If a tuple, mu_sigma must have 2 elements")
     elif mu_sigma is not None:
-        ncomp = algo_options.get('ncomp', None)
+        ncomp = algo_options.get('ncomp', 1)
         annulus_width = algo_options.get('annulus_width', int(fwhm))
         if weights is not None:
             if not len(weights) == cube.shape[0]:
@@ -280,7 +280,7 @@ def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo,
     # Calculate 1 sigma of distribution of deviations
     print(offset.shape)
     if force_rPA:
-        offset = offset[:, 2]
+        offset = offset[:, 2:]
         print(offset.shape)
     if trim_outliers:
         std = np.std(offset, axis=0)
@@ -293,12 +293,16 @@ def speckle_noise_uncertainty(cube, p_true, angle_range, derot_angles, algo,
     if bins is None:
         bins = int(offset.shape[0]/10)
 
+    if force_rPA:
+        labels = []
+    else:
+        labels = ['r', 'theta']
+
     if cube.ndim == 3:
-        labels=['r', 'theta', 'f']
+        labels.append('f')
     else:
-        labels=['r', 'theta']
         for ch in range(nch):
-            labels.append('f{}'.append(ch))
+            labels.append('f{}'.format(ch))
 
     mean_dev, sp_unc = confidence(offset, cfd=68.27, bins=bins, 
                                   gaussian_fit=True, verbose=True, save=False,
@@ -327,14 +331,14 @@ def _estimate_speckle_one_angle(angle, cube_pf, psfn, angs, r_true, f_true,
                                      imlib=imlib, interpolation=interpolation,
                                      verbose=False)
 
-    ncomp = algo_options.get('ncomp', None)
+    ncomp = algo_options.get('ncomp', 1)
     annulus_width = algo_options.get('annulus_width', int(fwhm))
 
     if cube_pf.ndim == 4:
         p_ini = [r_true, angle]
         for f in f_true:
             p_ini.append(f)
-            p_ini = tuple(p_ini)
+        p_ini = tuple(p_ini)
     else:
         p_ini = (r_true, angle, f_true)
 

vip_hci/fm/scattered_light_disk.py

@@ -587,6 +587,7 @@ def set_radial_density(self, ain=5., aout=-5., a=60.,
                 Gamma_out = self.aout+self.beta
                 self.apeak_surface_density = self.a * np.power(-Gamma_in/Gamma_out,
                                                                1./(2.*(Gamma_in-Gamma_out)))
+                # the above formula comes from Augereau et al. 1999.
             else:
                 self.apeak = self.a
                 self.apeak_surface_density = self.a
@@ -668,6 +669,11 @@ def half_max_density(r): return rad_density(r) / \
                     a_plus_hwhm,
                     a_plus_hwhm -
                     a_minus_hwhm))
+            msg4 = 'Semi-major axis at maximum dust surface density: {0:.1f}au ' \
+                   '(same as ref sma if ain=-aout)'
+            print(
+                msg4.format(
+                    self.apeak_surface_density))
             print('Ellipse p parameter: {0:.1f}au'.format(self.p))
         print('Ellipticity: {0:.3f}'.format(self.e))
         print('Inner slope: {0:.2f}'.format(self.ain))

vip_hci/invprob/fmmf.py

@@ -470,7 +470,7 @@ def _perturb(frame, model_matrix, numbasis, evals_matrix, evecs_matrix,
     """
     Function allowing the estimation of the PSF forward model when relying on
     KLIP for the computation of the speckle field. The code is based on the
-    PyKLIP library considering only the ADI case with a singlle number of
+    PyKLIP library considering only the ADI case with a single number of
     principal components considered. For more details about the code, consider
     the PyKLIP library or the original articles (Pueyo, L. 2016, ApJ, 824, 117
     or Ruffio, J.-B., Macintosh, B., Wang, J. J., & Pueyo, L. 2017, ApJ, 842)

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
+    return shifted_image

vip_hci/metrics/contrcurve.py

@@ -556,6 +556,7 @@ def throughput(cube, angle_list, psf_template, fwhm, algo, nbranch=1, theta=0,
     parangles = angle_list
     imlib = algo_dict.get('imlib', 'vip-fft')
     interpolation = algo_dict.get('interpolation', 'lanczos4')
+    scaling = algo_dict.get('scaling', None)
 
     if array.ndim != 3 and array.ndim != 4:
         raise TypeError('The input array is not a 3d or 4d cube')
@@ -645,9 +646,13 @@ def throughput(cube, angle_list, psf_template, fwhm, algo, nbranch=1, theta=0,
                                                       separation=fwhm_med,
                                                       fwhm=fwhm_med,
                                                       wedge=wedge)
-    noise_noscal, _, _ = noise_per_annulus(frame_nofc_noscal,
-                                           separation=fwhm_med, fwhm=fwhm_med,
-                                           wedge=wedge)
+    if scaling is not None:
+        noise_noscal, _, _ = noise_per_annulus(frame_nofc_noscal,
+                                               separation=fwhm_med, 
+                                               fwhm=fwhm_med, wedge=wedge)
+    else:
+        noise_noscal = noise.copy()
+
     vector_radd = vector_radd[inner_rad-1:]
     noise = noise[inner_rad-1:]
     res_level = res_level[inner_rad-1:]

vip_hci/preproc/badpixremoval.py

@@ -289,6 +289,7 @@ def cube_fix_badpix_isolated(array, bpm_mask=None, sigma_clip=3, num_neig=5,
             frame = array_out[i]
             smoothed = median_filter(frame, size, mode='mirror')
             frame[np.where(final_bpm)] = smoothed[np.where(final_bpm)]
+            array_out[i] = frame
             if verbose:
                 count_bp += np.sum(final_bpm)
 

vip_hci/preproc/cosmetics.py

@@ -139,7 +139,7 @@ def frame_pad(array, fac, fillwith=0, loc=0, scale=1, keep_parity=True,
 
     Parameters
     ----------
-    array : numpy ndarray
+    array : 2D numpy ndarray
         Input frame.
     fac : float > 1 or tuple of 2 floats > 1.
         Ratio of the size between padded and input frame. If a tuple, 
@@ -422,8 +422,8 @@ def approx_stellar_position(cube, fwhm, return_test=False, verbose=False):
     Returns
     -------
     star_approx_idx: 2d numpy array
-        Array of y and x approx coordinates of the star in each channel of the 
-        cube.
+        Array of y and x approximate coordinates of the star in each channel of 
+        the cube. Dimensions are nz x 2.
     test_result: 1d numpy array
         [return_test=True] It also returns the test result vector.
     """

vip_hci/preproc/derotation.py

@@ -134,6 +134,7 @@ def frame_rotate(array, angle, imlib='vip-fft', interpolation='lanczos4',
             _, med, stddev = sigma_clipped_stats(array_nan, sigma=1.5,
                                                  cenfunc=np.nanmedian,
                                                  stdfunc=np.nanstd)
+
         # pad and interpolate, about 1.2x original size
         if imlib == 'vip-fft':
             fac = 1.5
@@ -524,7 +525,7 @@ def _define_annuli(angle_list, ann, n_annuli, fwhm, radius_int, annulus_width,
 
 def rotate_fft(array, angle):
     """ Rotates a frame or 2D array using Fourier transform phases.
-    Rotation is equivalent to 3 consecutive linear shears, or 3 consecutive 
+    Rotation is equivalent to 3 consecutive linear shears, or 3 consecutive 1D
     FFT phase shifts. See details in [LAR97]_. 
 
     Parameters
@@ -590,7 +591,7 @@ def rotate_fft(array, angle):
     arr_y = arr_xy[0]-cy
     arr_x = arr_xy[1]-cx
 
-    # TODO: FFT padding not currently working properly. Only option '0' works.
+    # TODO: make FFT padding work for other option than '0'.
     s_x = _fft_shear(array_in, arr_x, a, ax=1, pad=0)
     s_xy = _fft_shear(s_x, arr_y, b, ax=0, pad=0)
     s_xyx = _fft_shear(s_xy, arr_x, a, ax=1, pad=0)

vip_hci/preproc/recentering.py

@@ -1330,7 +1330,7 @@ def cube_recenter_2dfit(array, xy=None, fwhm=4, subi_size=5, model='gauss',
         FWHM size in pixels, either one value (float) that will be the same for
         the whole cube, or an array of floats with the same dimension as the
         0th dim of array, containing the fwhm for each channel (e.g. in the case
-        of an ifs cube, where the fwhm varies with wavelength)
+        of an IFS cube, where the FWHM varies with wavelength).
     subi_size : int, optional
         Size of the square subimage sides in pixels.
     model : str, optional

vip_hci/psfsub/pca_fullfr.py

@@ -1022,8 +1022,8 @@ def _adi_rdi_pca(cube, cube_ref, angle_list, ncomp, scaling, mask_center_px,
     return pcs, recon, residuals_cube, residuals_cube_, frame
 
 
-def _project_subtract(cube, cube_ref, ncomp, scaling, mask_center_px,
-                      svd_mode, verbose, full_output, indices=None, frame=None,
+def _project_subtract(cube, cube_ref, ncomp, scaling, mask_center_px, svd_mode, 
+                      verbose, full_output, indices=None, frame=None,
                       cube_sig=None):
     """
     PCA projection and model PSF subtraction. Used as a helping function by

vip_hci/psfsub/pca_local.py

@@ -36,8 +36,8 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
                 delta_sep=(0.1, 1), ncomp=1, svd_mode='lapack', nproc=1,
                 min_frames_lib=2, max_frames_lib=200, tol=1e-1, scaling=None,
                 imlib='vip-fft', interpolation='lanczos4', collapse='median',
-                collapse_ifs='mean', ifs_collapse_range='all',
-                full_output=False, verbose=True, weights=None, cube_sig=None,
+                collapse_ifs='mean', ifs_collapse_range='all', theta_init=0, 
+                weights=None, cube_sig=None, full_output=False, verbose=True, 
                 **rot_options):
     """ PCA model PSF subtraction for ADI, ADI+RDI or ADI+mSDI (IFS) data. The
     PCA model is computed locally in each annulus (or annular sectors according
@@ -191,8 +191,10 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
         Whether to return the final median combined image only or with other
         intermediate arrays.
     verbose : bool, optional
-        If True prints to stdout intermediate info. If set to 2, the number of
-        frames in library and number of PCs are printed for each annular quadrant.
+        If True prints to stdout intermediate info.
+    theta_init : int
+        Initial azimuth [degrees] of the first segment, counting from the
+        positive x-axis counterclockwise (irrelevant if n_segments=1).
     weights: 1d numpy array or list, optional
         Weights to be applied for a weighted mean. Need to be provided if
         collapse mode is 'wmean'.
@@ -225,7 +227,7 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
                            n_segments, delta_rot, ncomp, svd_mode, nproc,
                            min_frames_lib, max_frames_lib, tol, scaling, imlib,
                            interpolation, collapse, True, verbose, cube_ref,
-                           weights, cube_sig, **rot_options)
+                           theta_init, weights, cube_sig, **rot_options)
 
         cube_out, cube_der, frame = res
         if full_output:
@@ -311,7 +313,8 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
         res = pool_map(nproc, _pca_sdi_fr, iterable(range(n)), scale_list,
                        radius_int, fwhm, asize, n_segments, delta_sep, ncomp,
                        svd_mode, tol, scaling, imlib, interpolation,
-                       collapse_ifs, ifs_collapse_range, verbose=verbose)
+                       collapse_ifs, ifs_collapse_range, theta_init, 
+                       verbose=verbose)
         residuals_cube_channels = np.array(res)
 
         # Exploiting rotational variability
@@ -338,8 +341,8 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
                                fwhm, asize, n_segments, delta_rot, ncomp2,
                                svd_mode, nproc, min_frames_lib, max_frames_lib,
                                tol, scaling, imlib, interpolation, collapse,
-                               full_output, verbose, None, weights, cube_sig,
-                               **rot_options)
+                               full_output, verbose, None, theta_init, weights, 
+                               cube_sig, **rot_options)
             if full_output:
                 cube_out, cube_der, frame = res
             else:
@@ -361,7 +364,7 @@ def pca_annular(cube, angle_list, cube_ref=None, scale_list=None, radius_int=0,
 
 def _pca_sdi_fr(fr, scal, radius_int, fwhm, asize, n_segments, delta_sep, ncomp,
                 svd_mode, tol, scaling, imlib, interpolation, collapse,
-                ifs_collapse_range):
+                ifs_collapse_range, theta_init):
     """ Optimized PCA subtraction on a multi-spectral frame (IFS data).
     """
     z, n, y_in, x_in = ARRAY.shape
@@ -402,7 +405,7 @@ def _pca_sdi_fr(fr, scal, radius_int, fwhm, asize, n_segments, delta_sep, ncomp,
         ann_center = inner_radius + (asize / 2)
 
         indices = get_annulus_segments(multispec_fr[0], inner_radius, asize,
-                                       n_segments[ann])
+                                       n_segments[ann], theta_init)
         # Library matrix is created for each segment and scaled if needed
         for seg in range(n_segments[ann]):
             yy = indices[seg][0]
@@ -442,8 +445,8 @@ def _pca_adi_rdi(cube, angle_list, radius_int=0, fwhm=4, asize=2, n_segments=1,
                  delta_rot=1, ncomp=1, svd_mode='lapack', nproc=None,
                  min_frames_lib=2, max_frames_lib=200, tol=1e-1, scaling=None,
                  imlib='vip-fft', interpolation='lanczos4', collapse='median',
-                 full_output=False, verbose=1, cube_ref=None, weights=None,
-                 cube_sig=None, **rot_options):
+                 full_output=False, verbose=1, cube_ref=None, theta_init=0, 
+                 weights=None, cube_sig=None, **rot_options):
     """ PCA exploiting angular variability (ADI fashion).
     """
     array = cube
@@ -501,7 +504,7 @@ def _pca_adi_rdi(cube, angle_list, radius_int=0, fwhm=4, asize=2, n_segments=1,
                                      n_segments_ann, verbose, True)
         pa_thr, inner_radius, ann_center = res_ann_par
         indices = get_annulus_segments(array[0], inner_radius, asize,
-                                       n_segments_ann)
+                                       n_segments_ann, theta_init)
         # Library matrix is created for each segment and scaled if needed
         for j in range(n_segments_ann):
             yy = indices[j][0]
@@ -588,6 +591,9 @@ def do_pca_patch(matrix, frame, angle_list, fwhm, pa_threshold, ann_center,
 
         if data_ref.shape[0] < min_frames_lib and matrix_ref is None:
             raise RuntimeError(msg.format(len(indices_left), min_frames_lib))
+    else:
+        data_ref = None
+
     if matrix_ref is not None:
         #data_ref = None
         # if matrix_ref is not None:
@@ -602,6 +608,13 @@ def do_pca_patch(matrix, frame, angle_list, fwhm, pa_threshold, ann_center,
         else:
             data_ref = matrix
 
+    if matrix_ref is not None and matrix_sig_segm is None:
+        # Stacking the ref and the target (pa thresh) libraries
+        if data_ref is not None:
+            data_ref = np.vstack((matrix_ref, data_ref))
+        else:
+            data_ref = matrix_ref
+
     curr_frame = matrix[frame]  # current frame
     if matrix_sig_segm is not None:
         curr_frame_emp = matrix[frame]-matrix_sig_segm[frame]

vip_hci/var/filters.py

@@ -136,7 +136,7 @@ def cube_filter_highpass(array, mode='laplacian', verbose=True, **kwargs):
     elif array.ndim == 4:
         for i in Progressbar(range(array.shape[1]), verbose=verbose):
             for lam in range(array.shape[0]):
-                array_out[lam][i] = frame_filter_highpass(array[lam][i],
+                array_out[lam][i] = frame_filter_highpass(array[lam][i], 
                                                           mode=mode, **kwargs)
     else:
         raise TypeError('Input array is not a 3d or 4d cube')