Skip to content

Commit

Permalink
Merge pull request #117 from carlgogo/master
Browse files Browse the repository at this point in the history 
VIP v0.8.8. Fixed multiprocessing for annular PCA, LLSG and ``snrmap`` functions in Python 3.6
  • Loading branch information
@carlos-gg
carlos-gg authored Feb 13, 2018
2 parents 7eec3f0 + 022bde0 commit c19f4ac
Show file tree
Hide file tree
Showing 21 changed files with 950 additions and 716 deletions.
20 changes: 5 additions & 15 deletions docs/source/faq.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,21 +23,11 @@ matplotlibrc file is located in:
$HOME/.config/matplotlib/matplotlibrc
* Why do I get, in OSX, the following RuntimeError?
*Python is not installed as a framework. The Mac OS X backend will not be able
to function correctly if Python is not installed as a framework. See the
Python documentation for more information on installing Python as a framework
on Mac OS X. Please either reinstall Python as a framework, or try one of the
other backends.*
Again, this is a matplotlib-backend related issue. Read the link in the previous
answer. It can be solved setting the backend to WXAgg or TkAgg. Optionally, you
can change your matplotlib backend **before** importing ``VIP``:

.. code-block:: python
import matplotlib as mpl
mpl.use('TkAgg')
import vip
* I get a RuntimeError (Python is not installed as a framework) in Macosx,
when using python3 and conda. Why?
This is caused by using matplotlib with virtual environments or conda in Macosx,
and has nothing to do with VIP. See 'Working with Matplotlib on OSX' in the
Matplotlib FAQ for more information.


* The ``VIP`` setup.py script doesn't finish the job, it seems to be stuck. What
Expand Down
2 changes: 1 addition & 1 deletion docs/source/readme.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ made by collaborators from several teams (take a look at the tab contributors on
it doesn't mean it's free from bugs. The code is continuously evolving and
therefore feedback/contributions are greatly appreciated. If you want to report
a bug, suggest or add a functionality please create an issue or send a pull
request on `this repository <https://github.com/vortex-exoplanet/VIP>`_.
request `here <https://github.com/vortex-exoplanet/VIP>`_.


Documentation
Expand Down
2 changes: 1 addition & 1 deletion readme.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ made by collaborators from several teams (take a look at the tab contributors on
it doesn't mean it's free from bugs. The code is continuously evolving and
therefore feedback/contributions are greatly appreciated. If you want to report
a bug, suggest or add a functionality please create an issue or send a pull
request on `this repository <https://github.com/vortex-exoplanet/VIP>`_.
request `here <https://github.com/vortex-exoplanet/VIP>`_.


Documentation
Expand Down
2 changes: 1 addition & 1 deletion vip_hci/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from . import stats
from . import var

__version__ = "0.8.7"
__version__ = "0.8.8"


print("---------------------------------------------------")
Expand Down
34 changes: 17 additions & 17 deletions vip_hci/llsg/llsg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ def llsg(cube, angle_list, fwhm, rank=10, thresh=1, max_iter=10,
residuals_tol=1e-1, cevr=0.9, thresh_mode='soft', nproc=1,
asize=None, n_segments=4, azimuth_overlap=None, radius_int=None,
random_seed=None, imlib='opencv', interpolation='lanczos4',
high_pass=False, collapse='median', full_output=False, verbose=True,
high_pass=None, collapse='median', full_output=False, verbose=True,
debug=False):
""" Local Low-rank plus Sparse plus Gaussian-noise decomposition (LLSG) as
described in Gomez Gonzalez et al. 2016. This first version of our algorithm
Expand Down Expand Up @@ -228,22 +228,22 @@ def llsg(cube, angle_list, fwhm, rank=10, thresh=1, max_iter=10,

else:
pool = Pool(processes=int(nproc))
res = pool.map(EFT, itt.izip(itt.repeat(_decompose_patch),
itt.repeat(indices),
range(n_segments_ann),
itt.repeat(n_segments_ann),
itt.repeat(rank),
itt.repeat(low_rank_ref),
itt.repeat(low_rank_mode),
itt.repeat(thresh),
itt.repeat(thresh_mode),
itt.repeat(max_iter),
itt.repeat(auto_rank_mode),
itt.repeat(cevr),
itt.repeat(residuals_tol),
itt.repeat(random_seed),
itt.repeat(debug),
itt.repeat(full_output)))
res = pool.map(EFT, zip(itt.repeat(_decompose_patch),
itt.repeat(indices),
range(n_segments_ann),
itt.repeat(n_segments_ann),
itt.repeat(rank),
itt.repeat(low_rank_ref),
itt.repeat(low_rank_mode),
itt.repeat(thresh),
itt.repeat(thresh_mode),
itt.repeat(max_iter),
itt.repeat(auto_rank_mode),
itt.repeat(cevr),
itt.repeat(residuals_tol),
itt.repeat(random_seed),
itt.repeat(debug),
itt.repeat(full_output)))
patches = np.array(res)
pool.close()

Expand Down
4 changes: 1 addition & 3 deletions vip_hci/negfc/mcmc_sampling.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,7 @@
from __future__ import print_function

__author__ = 'O. Wertz, Carlos Alberto Gomez Gonzalez'
__all__ = ['lnprior',
'lnlike',
'mcmc_negfc_sampling',
__all__ = ['mcmc_negfc_sampling',
'chain_zero_truncated',
'show_corner_plot',
'show_walk_plot',
Expand Down
4 changes: 3 additions & 1 deletion vip_hci/negfc/simplex_fmerit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,13 @@
"""
from __future__ import print_function

__all__ = []

import numpy as np
from skimage.draw import circle
from ..phot import cube_inject_companions
from ..var import frame_center
from ..pca import pca_annulus
from ..pca.utils_pca import pca_annulus


def chisquare(modelParameters, cube, angs, plsc, psfs_norm, fwhm, annulus_width,
Expand Down
4 changes: 1 addition & 3 deletions vip_hci/negfc/simplex_optim.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,7 @@
from ..var import frame_center
from ..conf import time_ini, timing, sep

__all__ = ['firstguess_from_coord',
'firstguess_simplex',
'firstguess']
__all__ = ['firstguess']



Expand Down
30 changes: 15 additions & 15 deletions vip_hci/pca/pca_fullfr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,8 @@ def pca(cube, angle_list=None, cube_ref=None, scale_list=None, ncomp=1, ncomp2=1
ncomp is the number of PCs from the library of spectral channels.
ncomp2 : int, optional
How many PCs are used for IFS+ADI datacubes in the second stage PCA.
ncomp2 goes up to the number of multi-spectral frames.
ncomp2 goes up to the number ADI frames. If None then the second stage
of the PCA is ignored and the residuals are derotated and combined.
svd_mode : {'lapack', 'arpack', 'eigen', 'randsvd', 'cupy', 'eigencupy', 'randcupy'}, str
Switch for the SVD method/library to be used. ``lapack`` uses the LAPACK
linear algebra library through Numpy and it is the most conventional way
Expand Down Expand Up @@ -193,7 +194,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
else:
ref_lib = matrix

if indices is not None and frame is not None: # one row (frame) at a time
if indices is not None and frame is not None: # one row (frame) at a time
ref_lib = ref_lib[indices]
if ref_lib.shape[0] <= 10:
msg = 'Too few frames left in the PCA library (<10). '
Expand All @@ -209,7 +210,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
return ref_lib.shape[0], residuals, reconstructed
else:
return ref_lib.shape[0], residuals
else: # the whole matrix
else: # the whole matrix
V = svd_wrapper(ref_lib, svd_mode, ncomp, debug, verbose)

if verbose: timing(start_time)
Expand Down Expand Up @@ -263,7 +264,8 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
if not scale_list.shape[0]==cube.shape[0]:
raise TypeError('Scaling factors vector has wrong length')

if verbose: start_time = time_ini()
if verbose:
start_time = time_ini()

if check_mem:
input_bytes = cube.nbytes
Expand Down Expand Up @@ -322,7 +324,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
# SDI (IFS) + ADI: cube with multiple spectral channels + rotation
# shape of cube: [# channels, # adi-frames, Y, X]
#***********************************************************************
elif cube.ndim==4 and angle_list is not None:
elif cube.ndim == 4 and angle_list is not None:
if verbose:
print('{:} spectral channels in IFS cube'.format(z))
print('{:} ADI frames in all channels'.format(n))
Expand Down Expand Up @@ -353,22 +355,20 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
residuals_cube_channels[i] = frame
bar.update()

# de-rotation of the PCA processed channels
if ncomp2 > z:
ncomp2 = min(ncomp2, z)
msg = 'Number of PCs too high (max PCs={}), using instead {:} PCs.'
print(msg.format(n, ncomp2))

elif ncomp2 is None:
# de-rotation of the PCA processed channels, ADI fashion
if ncomp2 is None:
residuals_cube_channels_ = cube_derotate(residuals_cube_channels,
angle_list, imlib=imlib,
interpolation=interpolation)
angle_list, imlib=imlib,
interpolation=interpolation)
frame = cube_collapse(residuals_cube_channels_, mode=collapse)
if verbose:
print('De-rotating and combining')
timing(start_time)

else:
if ncomp2 > n:
ncomp2 = min(ncomp2, n)
msg = 'Number of PCs too high (max PCs={}), using instead {:} PCs.'
print(msg.format(n, ncomp2))
res_ifs_adi = subtract_projection(residuals_cube_channels, None,
ncomp2, scaling, mask_center_px,
debug, svd_mode, False,
Expand Down
68 changes: 28 additions & 40 deletions vip_hci/pca/pca_local.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,11 @@
from ..stats import descriptive_stats



def pca_rdi_annular(cube, angle_list, cube_ref, radius_int=0, asize=1,
ncomp=1, svd_mode='randsvd', min_corr=0.9, fwhm=4,
scaling='temp-standard', imlib='opencv',
interpolation='lanczos4', collapse='median',
full_output=False, verbose=True, debug=False):
full_output=False, verbose=True):
""" Annular PCA with Reference Library + Correlation + standardization
In the case of having a large number of reference images, e.g. for a survey
Expand Down Expand Up @@ -87,9 +86,7 @@ def pca_rdi_annular(cube, angle_list, cube_ref, radius_int=0, asize=1,
Whether to return the final median combined image only or with other
intermediate arrays.
verbose : {True, False}, bool optional
If True prints to stdout intermediate info.
debug : {False, True}, bool optional
Whether to output some intermediate information.
If True prints to stdout intermediate info.
Returns
-------
Expand Down Expand Up @@ -204,7 +201,7 @@ def do_pca_annulus(ncomp, matrix, svd_mode, noise_error, data_ref):
return frame


def pca_adi_annular(cube, angle_list, radius_int=0, fwhm=4, asize=3,
def pca_adi_annular(cube, angle_list, radius_int=0, fwhm=4, asize=3,
delta_rot=1, ncomp=1, svd_mode='randsvd', nproc=1,
min_frames_pca=10, tol=1e-1, scaling=None, quad=False,
imlib='opencv', interpolation='lanczos4', collapse='median',
Expand All @@ -219,10 +216,7 @@ def pca_adi_annular(cube, angle_list, radius_int=0, fwhm=4, asize=3,
algebra calcularions, which comes by default in every OSX system, is broken
for multiprocessing. Avoid using this function unless you have compiled
Python against other linear algebra library. An easy fix is to install
latest ANACONDA (2.5 or later) distribution which ships MKL library
(replacing the problematic ACCELERATE). On linux with the default
LAPACK/BLAS libraries it successfully distributes the processes among all
the existing cores.
(ana)conda and the openblas or MKL libraries (replacing ACCELERATE).
Parameters
----------
Expand Down Expand Up @@ -378,10 +372,10 @@ def pca_adi_annular(cube, angle_list, radius_int=0, fwhm=4, asize=3,
#***************************************************************
# We loop the frames and do the PCA to obtain the residuals cube
#***************************************************************
residuals = do_pca_loop(matrix_quad, yy, xx, nproc, angle_list,
fwhm, pa_thr, scaling, ann_center,
svd_mode, ncompann, min_frames_pca, tol,
debug, verbose)
residuals = do_pca_loop(matrix_quad, nproc, angle_list, fwhm,
pa_thr, scaling, ann_center, svd_mode,
ncompann, min_frames_pca, tol, debug,
verbose)

for frame in range(n):
cube_out[frame][yy, xx] = residuals[frame]
Expand All @@ -404,9 +398,9 @@ def pca_adi_annular(cube, angle_list, radius_int=0, fwhm=4, asize=3,
#*******************************************************************
# We loop the frames and do the PCA to obtain the residuals cube
#*******************************************************************
residuals = do_pca_loop(matrix_ann, yy, xx, nproc, angle_list, fwhm,
pa_thr, scaling, ann_center, svd_mode,
ncompann, min_frames_pca, tol, debug, verbose)
residuals = do_pca_loop(matrix_ann, nproc, angle_list, fwhm, pa_thr,
scaling, ann_center, svd_mode, ncompann,
min_frames_pca, tol, debug, verbose)

for frame in range(n):
cube_out[frame][yy, xx] = residuals[frame]
Expand Down Expand Up @@ -504,11 +498,11 @@ def find_indices(angle_list, frame, thr, truncate):
else:
half2 = range(index_foll, min(n, int(thr/2+index_foll)))
half1 = range(max(0,index_prev-thr+len(half2)), index_prev)
return np.array(half1+half2)
return np.array(list(half1) + list(half2))



def do_pca_loop(matrix, yy, xx, nproc, angle_list, fwhm, pa_threshold, scaling,
def do_pca_loop(matrix, nproc, angle_list, fwhm, pa_threshold, scaling,
ann_center, svd_mode, ncomp, min_frames_pca, tol, debug, verbose):
"""
"""
Expand All @@ -521,40 +515,34 @@ def do_pca_loop(matrix, yy, xx, nproc, angle_list, fwhm, pa_threshold, scaling,
#***************************************************************
ncomps = []
nfrslib = []
if nproc==1:
if nproc == 1:
residualarr = []
for frame in range(n):
res = do_pca_patch(matrix_ann, frame, angle_list, fwhm,
pa_threshold, scaling, ann_center,
svd_mode, ncomp, min_frames_pca, tol,
debug)
pa_threshold, ann_center, svd_mode, ncomp,
min_frames_pca, tol, debug)
residualarr.append(res[0])
ncomps.append(res[1])
nfrslib.append(res[2])
residuals = np.array(residualarr)

elif nproc>1:
elif nproc > 1:
#***********************************************************
# A multiprocessing pool is created to process the frames in
# a parallel way. SVD/PCA is done in do_pca_patch function
#***********************************************************
pool = Pool(processes=int(nproc))
res = pool.map(EFT, itt.izip(itt.repeat(do_pca_patch),
itt.repeat(matrix_ann),
range(n), itt.repeat(angle_list),
itt.repeat(fwhm),
itt.repeat(pa_threshold),
itt.repeat(scaling),
itt.repeat(ann_center),
itt.repeat(svd_mode),
itt.repeat(ncomp),
itt.repeat(min_frames_pca),
itt.repeat(tol),
itt.repeat(debug)))
res = pool.map(EFT, zip(itt.repeat(do_pca_patch),
itt.repeat(matrix_ann), range(n),
itt.repeat(angle_list), itt.repeat(fwhm),
itt.repeat(pa_threshold),
itt.repeat(ann_center), itt.repeat(svd_mode),
itt.repeat(ncomp), itt.repeat(min_frames_pca),
itt.repeat(tol)))
res = np.array(res)
residuals = np.array(res[:,0])
ncomps = res[:,1]
nfrslib = res[:,2]
ncomps = res[:, 1]
nfrslib = res[:, 2]
pool.close()

# number of frames in library printed for each annular quadrant
Expand All @@ -567,8 +555,8 @@ def do_pca_loop(matrix, yy, xx, nproc, angle_list, fwhm, pa_threshold, scaling,
return residuals


def do_pca_patch(matrix, frame, angle_list, fwhm, pa_threshold, scaling,
ann_center, svd_mode, ncomp, min_frames_pca, tol, debug):
def do_pca_patch(matrix, frame, angle_list, fwhm, pa_threshold, ann_center,
svd_mode, ncomp, min_frames_pca, tol):
"""
Does the SVD/PCA for each frame patch (small matrix). For each frame we
find the frames to be rejected depending on the amount of rotation. The
Expand Down
Loading

0 comments on commit c19f4ac

Please sign in to comment.