Merge pull request #583 from Sand-jrd/master

Add "Temporal" option to pca fullfr

tests/test_pipeline_adi.py

@@ -77,6 +77,9 @@ def algo_nmf_annular(ds):
 def algo_pca(ds):
     return vip.psfsub.pca(ds.cube, ds.angles, svd_mode='arpack')
 
+def algo_pca_left_eigv(ds):
+    return vip.psfsub.pca(ds.cube, ds.angles, left_eigv=True)
+
 def algo_pca_linalg(ds):
     return vip.psfsub.pca(ds.cube, ds.angles, svd_mode='eigen')
 
@@ -163,6 +166,7 @@ def verify_expcoord(vectory, vectorx, exp_yx):
                  (algo_frdiff, snrmap_fast),
                  (algo_frdiff4, snrmap_fast),
                  (algo_pca, snrmap_fast),
+                 (algo_pca_left_eigv, snrmap_fast),
                  (algo_pca_linalg, snrmap_fast),
                  (algo_pca_drot, snrmap_fast),
                  (algo_pca_cevr, snrmap_fast),

tests/test_pipeline_sdi.py

@@ -94,11 +94,17 @@ def algo_pca_single(ds, sc):
     return vip.psfsub.pca(ds.cube, ds.angles, scale_list=sc,
                           adimsdi='single', ncomp=10)
 
+def algo_pca_single_left_eigv(ds, sc):
+    return vip.psfsub.pca(ds.cube, ds.angles, scale_list=sc,
+                          adimsdi='single', ncomp=10, left_eigv=True)
 
 def algo_pca_double(ds, sc):
     return vip.psfsub.pca(ds.cube, ds.angles, scale_list=sc,
                           adimsdi='double', ncomp=(1, 2))
 
+def algo_pca_double_left_eigv(ds, sc):
+    return vip.psfsub.pca(ds.cube, ds.angles, scale_list=sc,
+                          adimsdi='double', ncomp=(1, 2), left_eigv=True)
 
 def algo_pca_annular(ds, sc):
     return vip.psfsub.pca_annular(ds.cube, ds.angles, scale_list=sc,
@@ -156,7 +162,9 @@ def verify_expcoord(vectory, vectorx, exp_yx):
                  (algo_xloci, snrmap_fast),
                  (algo_xloci_double, snrmap_fast),
                  (algo_pca_single, snrmap_fast),
+                 (algo_pca_single_left_eigv, snrmap_fast),
                  (algo_pca_double, snrmap_fast),
+                 (algo_pca_double_left_eigv, snrmap_fast),
                  (algo_pca_annular, None),
                  ],
              ids=lambda x: (x.__name__.replace("algo_", "") if callable(x) else x))

vip_hci/psfsub/pca_fullfr.py

@@ -97,6 +97,7 @@ def pca(
     verbose=True,
     weights=None,
     conv=False,
+    left_eigv=False,
     cube_sig=None,
     **rot_options
 ):
@@ -312,6 +313,9 @@ def pca(
     weights: 1d numpy array or list, optional
         Weights to be applied for a weighted mean. Need to be provided if
         collapse mode is 'wmean'.
+    left_eigv : bool, optional
+        Whether to use rather left or right singularvectors
+        This mode is not compatible with 'mask_rdi' and 'batch'
     cube_sig: numpy ndarray, opt
         Cube with estimate of significant authentic signals. If provided, this
         will subtracted before projecting cube onto reference cube.
@@ -320,7 +324,7 @@ def pca(
         "edge_blend", "interp_zeros", "ker" (see documentation of
         ``vip_hci.preproc.frame_rotate``)
 
-    Returns
+    Return
     -------
     frame : numpy ndarray
         2D array, median combination of the de-rotated/re-scaled residuals cube.
@@ -371,6 +375,12 @@ def pca(
                 "with the full path on disk"
             )
 
+    if left_eigv : 
+        if (batch is not None or mask_rdi is not None or cube_ref is not None):
+            raise NotImplementedError( "left_eigv is not compatible"
+                                      "with 'mask_rdi' nor 'batch'"
+            )
+
     # checking memory (if in-memory numpy array is provided)
     if not isinstance(cube, str):
         input_bytes = cube_ref.nbytes if cube_ref is not None else cube.nbytes
@@ -416,6 +426,7 @@ def pca(
                 conv,
                 mask_rdi,
                 cube_sig,
+                left_eigv,
                 **rot_options
             )
             residuals_cube_channels, residuals_cube_channels_, frame = res_pca
@@ -443,6 +454,7 @@ def pca(
                 batch,
                 full_output=True,
                 weights=weights,
+                left_eigv=left_eigv,
                 **rot_options
             )
             if isinstance(ncomp, (int, float)):
@@ -534,6 +546,7 @@ def pca(
                     True,
                     weights,
                     cube_sig,
+                    left_eigv,
                     **rot_options
                 )
                 if batch is None:
@@ -630,6 +643,7 @@ def pca(
             True,
             weights,
             cube_sig,
+            left_eigv,
             **rot_options
         )
 
@@ -746,6 +760,7 @@ def _adi_pca(
     full_output,
     weights=None,
     cube_sig=None,
+    left_eigv=False,
     **rot_options
 ):
     """Handle the ADI PCA post-processing."""
@@ -803,14 +818,15 @@ def _adi_pca(
                     verbose,
                     full_output,
                     cube_sig=cube_sig,
+                    left_eigv=left_eigv,
                 )
                 if verbose:
                     timing(start_time)
                 if full_output:
                     residuals_cube = residuals_result[0]
                     reconstructed = residuals_result[1]
                     V = residuals_result[2]
-                    pcs = reshape_matrix(V, y, x)
+                    pcs = reshape_matrix(V, y, x) if not left_eigv else V.T
                     recon = reshape_matrix(reconstructed, y, x)
                 else:
                     residuals_cube = residuals_result
@@ -857,6 +873,7 @@ def _adi_pca(
                         ind,
                         frame,
                         cube_sig=cube_sig,
+                        left_eigv=left_eigv
                     )
                     if full_output:
                         nfrslib.append(res_result[0])
@@ -944,6 +961,7 @@ def _adimsdi_singlepca(
     batch,
     full_output,
     weights=None,
+    left_eigv=False,
     **rot_options
 ):
     """Handle the full-frame ADI+mSDI single PCA post-processing."""
@@ -1004,8 +1022,8 @@ def _adimsdi_singlepca(
         # When ncomp is a int/float and batch is None, standard ADI-PCA is run
         else:
             res_cube = _project_subtract(
-                big_cube, None, ncomp, scaling, mask_center_px, svd_mode, verbose, False
-            )
+                big_cube, None, ncomp, scaling, mask_center_px, svd_mode, 
+                verbose, False, left_eigv=left_eigv)
 
         if verbose:
             timing(start_time)
@@ -1112,6 +1130,7 @@ def _adimsdi_doublepca(
     conv=False,
     mask_rdi=None,
     cube_sig=None,
+    left_eigv=False,
     **rot_options
 ):
     """Handle the full-frame ADI+mSDI double PCA post-processing."""
@@ -1170,6 +1189,7 @@ def _adimsdi_doublepca(
         fwhm,
         conv,
         mask_rdi,
+        left_eigv,
     )
     residuals_cube_channels = np.array(res)
 
@@ -1212,6 +1232,7 @@ def _adimsdi_doublepca(
             verbose=False,
             full_output=False,
             cube_sig=cube_sig,
+            left_eigv=left_eigv,
         )
         if verbose:
             print("De-rotating and combining residuals")
@@ -1244,6 +1265,8 @@ def _adimsdi_doublepca_ifs(
     fwhm,
     conv,
     mask_rdi=None,
+    left_eigv=False,
+
 ):
     """Call by _adimsdi_doublepca with pool_map."""
     global ARRAY
@@ -1280,6 +1303,7 @@ def _adimsdi_doublepca_ifs(
                 svd_mode,
                 verbose=False,
                 full_output=False,
+                left_eigv=left_eigv,
             )
         else:
             residuals = np.zeros_like(cube_resc)
@@ -1402,6 +1426,7 @@ def _project_subtract(
     indices=None,
     frame=None,
     cube_sig=None,
+    left_eigv=False,
 ):
     """
     PCA projection and model PSF subtraction.
@@ -1426,6 +1451,8 @@ def _project_subtract(
         Verbosity.
     full_output : bool
         Whether to return intermediate arrays or not.
+    left_eigv : bool, optional
+        Whether to use rather left or right singularvectors
     indices : list
         Indices to be used to discard frames (a rotation threshold is used).
     frame : int
@@ -1446,8 +1473,9 @@ def _project_subtract(
     reconstructed : numpy ndarray
         [full_output=True] The reconstructed array.
     V : numpy ndarray
-        [full_output=True, indices is None, frame is None] The right singular
-        vectors of the input matrix, as returned by ``svd/svd_wrapper()``
+        [full_output=True, indices is None, frame is None]
+        The right singular vectors of the input matrix, as returned by 
+        ``svd/svd_wrapper()``
     """
     _, y, x = cube.shape
     if isinstance(ncomp, (int, np.int_)):
@@ -1474,6 +1502,10 @@ def _project_subtract(
             ref_lib = prepare_matrix(
                 cube_ref, scaling, mask_center_px, mode="fullfr", verbose=verbose
             )
+        elif left_eigv:
+            ref_lib = prepare_matrix(cube, scaling, mask_center_px, 
+                      mode="fullfr", verbose=verbose, discard_mask_pix=True
+            )
         else:
             ref_lib = matrix_emp
 
@@ -1487,22 +1519,36 @@ def _project_subtract(
                 )
             curr_frame = matrix[frame]  # current frame
             curr_frame_emp = matrix_emp[frame]
-            V = svd_wrapper(ref_lib, svd_mode, ncomp, False)
-            transformed = np.dot(curr_frame_emp, V.T)
-            reconstructed = np.dot(transformed.T, V)
+            if left_eigv :
+                V = svd_wrapper(ref_lib, svd_mode, ncomp, False, left_eigv=left_eigv)
+                transformed = np.dot(curr_frame_emp.T, V)
+                reconstructed = np.dot(V, transformed.T)
+            else :
+                V = svd_wrapper(ref_lib, svd_mode, ncomp, False)
+                transformed = np.dot(curr_frame_emp, V.T)
+                reconstructed = np.dot(transformed.T, V)
+
             residuals = curr_frame - reconstructed
+
             if full_output:
                 return ref_lib.shape[0], residuals, reconstructed
             else:
                 return ref_lib.shape[0], residuals
 
         # the whole matrix is processed at once
         else:
-            V = svd_wrapper(ref_lib, svd_mode, ncomp, verbose)
-            transformed = np.dot(V, matrix_emp.T)
-            reconstructed = np.dot(transformed.T, V)
+            if left_eigv :
+                V = svd_wrapper(ref_lib, svd_mode, ncomp, verbose, left_eigv=left_eigv)
+                transformed = np.dot(matrix_emp.T, V)
+                reconstructed = np.dot(V, transformed.T)
+            else :
+                V = svd_wrapper(ref_lib, svd_mode, ncomp, verbose)
+                transformed = np.dot(V, matrix_emp.T)
+                reconstructed = np.dot(transformed.T, V)
+
             residuals = matrix - reconstructed
             residuals_res = reshape_matrix(residuals, y, x)
+
             if full_output:
                 return residuals_res, reconstructed, V
             else:

vip_hci/psfsub/svd.py

@@ -335,7 +335,7 @@ def cevr_to_ncomp(self, cevr=0.9):
 
 
 def svd_wrapper(matrix, mode, ncomp, verbose, full_output=False,
-                random_state=None, to_numpy=True):
+                random_state=None, to_numpy=True, left_eigv=False):
     """ Wrapper for different SVD libraries (CPU and GPU).
 
     Parameters
@@ -384,16 +384,19 @@ def svd_wrapper(matrix, mode, ncomp, verbose, full_output=False,
     verbose: bool
         If True intermediate information is printed out.
     full_output : bool optional
-        If True the 3 terms of the SVD factorization are returned. If ``mode``
-        is eigen then only S and V are returned.
+        If True the 3 terms of the SVD factorization are returned. 
     random_state : int, RandomState instance or None, optional
         If int, random_state is the seed used by the random number generator.
         If RandomState instance, random_state is the random number generator.
         If None, the random number generator is the RandomState instance used
         by np.random. Used for ``randsvd`` mode.
     to_numpy : bool, optional
+        Whether to return intermediate arrays or not.
         If True (by default) the arrays computed in GPU are transferred from
         VRAM and converted to numpy ndarrays.
+    left_eigv : bool, optional
+        Whether to use rather left or right singularvectors
+        
 
     Returns
     -------
@@ -449,6 +452,9 @@ def svd_wrapper(matrix, mode, ncomp, verbose, full_output=False,
         for i in range(V.shape[1]):
             V[:, i] /= S    # scaling EVs by the square root of EVals
         V = V[:ncomp]
+        if full_output or left_eigv:
+            U = EV/np.sqrt(np.abs(e)) 
+            U = U[:ncomp]
         if verbose:
             print('Done PCA with numpy linalg eigh functions')
 
@@ -521,6 +527,10 @@ def svd_wrapper(matrix, mode, ncomp, verbose, full_output=False,
         V = V[:ncomp]
         if to_numpy:
             V = cupy.asnumpy(V)
+        if full_output or left_eigv:
+            U = EV/np.sqrt(np.abs(e)) 
+            U = U[:ncomp]
+            if to_numpy: U = cupy.asnumpy(U)
         if verbose:
             print('Done PCA with cupy eigh function (GPU)')
 
@@ -552,6 +562,10 @@ def svd_wrapper(matrix, mode, ncomp, verbose, full_output=False,
         V = V[:ncomp]
         if to_numpy:
             V = np.array(V)
+        if full_output or left_eigv:
+            U = EV/np.sqrt(np.abs(e)) 
+            U = U[:ncomp]
+            if to_numpy: U = cupy.asnumpy(U)    
         if verbose:
             print('Done PCA with pytorch eig function')
 
@@ -577,10 +591,15 @@ def svd_wrapper(matrix, mode, ncomp, verbose, full_output=False,
                 return V.T, S, U.T
             else:
                 return torch.transpose(V, 0, 1), S, torch.transpose(U, 0, 1)
-        elif mode in ('eigen', 'eigencupy', 'eigenpytorch'):
-            return S, V
         else:
             return U, S, V
+    elif left_eigv:
+        if mode == 'lapack':
+            return V.T
+        elif mode == 'pytorch':
+            return V
+        else:
+            return U
     else:
         if mode == 'lapack':
             return U.T