Skip to content

Commit

Permalink
Merge pull request #68 from carlgogo/master
Browse files Browse the repository at this point in the history 
VIP v0.6.2
  • Loading branch information
@carlos-gg
carlos-gg authored Dec 15, 2016
2 parents 60d58d5 + b062766 commit b1db7c6
Show file tree
Hide file tree
Showing 6 changed files with 129 additions and 135 deletions.
2 changes: 1 addition & 1 deletion vip/__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.6.1"
__version__ = "0.6.2"

print("---------------------------------------------------")
print(" oooooo oooo ooooo ooooooooo. ")
Expand Down
2 changes: 1 addition & 1 deletion vip/conf/mem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ def check_enough_memory(input_bytes, factor=1, verbose=True):
----------
input_bytes : float
The size in bytes of the inputs of a given function.
factor : int
factor : float
Scales how much memory is needed in terms of the size of input_bytes.
"""
mem = virtual_memory()
Expand Down
96 changes: 47 additions & 49 deletions vip/pca/pca_fullfr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,11 +212,11 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
if not cube.ndim>2:
raise TypeError('Input array is not a 3d or 4d array')
if angle_list is not None:
if (cube.ndim==3 and not (cube.shape[0] == angle_list.shape[0])) or \
if (cube.ndim==3 and not (cube.shape[0] == angle_list.shape[0])) or \
(cube.ndim==4 and not (cube.shape[1] == angle_list.shape[0])):
msg = "Angle list vector has wrong length. It must equal the number"
msg += " frames in the cube."
raise TypeError(msg)
msg = "Angle list vector has wrong length. It must equal the number"
msg += " frames in the cube."
raise TypeError(msg)
if source_xy is not None and delta_rot is None or fwhm is None:
msg = 'Delta_rot or fwhm parameters missing. They are needed for the '
msg += 'PA-based rejection of frames from the library'
Expand Down Expand Up @@ -268,15 +268,15 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
if ncomp > z:
ncomp = min(10, z)
msg = 'Number of PCs too high (max PCs={}), using instead {:} PCs.'
print msg.format(z, ncomp)
print(msg.format(z, ncomp))
scale_list = check_scal_vector(scale_list)

#***********************************************************************
# RDI (IFS): case of 3d cube with multiple spectral channels
#***********************************************************************
if cube.ndim==3:
if verbose:
print '{:} spectral channels in IFS cube'.format(z)
print('{:} spectral channels in IFS cube'.format(z))
# cube has been re-scaled to have the planets moving radially
cube, _, y, x, _, _ = scale_cube_for_pca(cube, scale_list)
residuals_result = subtract_projection(cube, None, ncomp, scaling,
Expand All @@ -299,7 +299,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
full_output=full_output,
inverse=True, y_in=y_in, x_in=x_in)
if verbose:
print 'Done re-scaling and combining'
print('Done re-scaling and combining')
timing(start_time)

#***********************************************************************
Expand All @@ -308,13 +308,13 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
#***********************************************************************
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)
print('{:} spectral channels in IFS cube'.format(z))
print('{:} ADI frames in all channels'.format(n))
residuals_cube_channels = np.zeros((n, y_in, x_in))

bar = pyprind.ProgBar(n, stream=1,
title='Looping through ADI frames')
for i in xrange(n):
for i in range(n):
cube_res, _, y, x, _, _ = scale_cube_for_pca(cube[:,i,:,:],
scale_list)
residuals_result = subtract_projection(cube_res, None, ncomp,
Expand All @@ -341,15 +341,14 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
if ncomp2 > n:
ncomp2 = min(10, n)
msg = 'Number of PCs too high (max PCs={}), using instead {:} PCs.'
print msg.format(n, ncomp2)
print(msg.format(n, ncomp2))

elif ncomp2 is None:
residuals_cube_channels_ = cube_derotate(residuals_cube_channels,
angle_list)
frame = cube_collapse(residuals_cube_channels_, mode=collapse)
if verbose:
msg = 'De-rotating and combining'
print msg
print('De-rotating and combining')
timing(start_time)

else:
Expand All @@ -360,9 +359,9 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
residuals_cube_channels_ = cube_derotate(res_ifs_adi, angle_list)
frame = cube_collapse(residuals_cube_channels_, mode=collapse)
if verbose:
msg = 'Done PCA per ADI multi-spectral frame, de-rotating and '
msg += 'combining'
print msg
msg = 'Done PCA per ADI multi-spectral frame, de-rotating '
msg += 'and combining'
print(msg)
timing(start_time)

#***************************************************************************
Expand All @@ -372,7 +371,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
if ncomp > n:
ncomp = min(ncomp,n)
msg = 'Number of PCs too high (max PCs={}), using instead {:} PCs.'
print msg.format(n, ncomp)
print(msg.format(n, ncomp))
residuals_result = subtract_projection(cube, cube_ref, ncomp, scaling,
mask_center_px, debug, svd_mode,
verbose, full_output)
Expand All @@ -388,7 +387,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
frame = cube_collapse(residuals_cube_, mode=collapse)

if verbose:
print 'Done de-rotating and combining'
print('Done de-rotating and combining')
timing(start_time)

#***************************************************************************
Expand All @@ -398,7 +397,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
if ncomp > n:
ncomp = min(ncomp,n)
msg = 'Number of PCs too high (max PCs={}), using instead {:} PCs.'
print msg.format(n, ncomp)
print(msg.format(n, ncomp))

if source_xy is None:
residuals_result = subtract_projection(cube, None, ncomp, scaling,
Expand All @@ -425,10 +424,10 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
new_pa_th = float(mid_range - mid_range * 0.1)
if verbose:
msg = 'PA threshold {:.2f} is too big, will be set to {:.2f}'
print msg.format(pa_thr, new_pa_th)
print(msg.format(pa_thr, new_pa_th))
pa_thr = new_pa_th

for frame in xrange(n):
for frame in range(n):
if ann_center > fwhm*3: # TODO: 3 optimal value? new parameter?
ind = find_indices(angle_list, frame, pa_thr, True)
else:
Expand Down Expand Up @@ -457,7 +456,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px,
residuals_cube_ = cube_derotate(residuals_cube, angle_list)
frame = cube_collapse(residuals_cube_, mode=collapse)
if verbose:
print 'Done de-rotating and combining'
print('Done de-rotating and combining')
timing(start_time)

if full_output and cube.ndim<4:
Expand Down Expand Up @@ -628,8 +627,8 @@ def grid(matrix, angle_list, y, x, mode, V, fwhm, fmerit, step, inti, intf,
if full_output: frlist = []
counter = 0
if debug:
print 'Step current grid:', step
print 'PCs | SNR'
print('Step current grid:', step)
print('PCs | SNR')
for pc in range(inti, intf+1, step):
if full_output:
snr, flux, frame = get_snr(matrix, angle_list, y, x, mode, V,
Expand All @@ -645,22 +644,23 @@ def grid(matrix, angle_list, y, x, mode, V, fwhm, fmerit, step, inti, intf,
if full_output: frlist.append(frame)
nsteps += 1
if truncate and nsteps>2 and snr<min_snr:
if debug: print 'SNR too small'
if debug: print('SNR too small')
break
if debug: print '{} {:.3f}'.format(pc, snr)
if debug: print('{} {:.3f}'.format(pc, snr))
if truncate and counter==5: break
argm = np.argmax(snrlist)

if len(pclist)==2: pclist.append(pclist[-1]+1)

if debug:
print 'Finished current stage'
print('Finished current stage')
try:
pclist[argm+1]
print 'Interval for next grid: ', pclist[argm-1], 'to', pclist[argm+1]
print('Interval for next grid: ', pclist[argm-1], 'to',
pclist[argm+1])
except:
print 'The optimal SNR seems to be outside of the given PC range'
print
print('The optimal SNR seems to be outside of the given PC range')
print()

if argm==0: argm = 1
if full_output:
Expand Down Expand Up @@ -731,10 +731,10 @@ def grid(matrix, angle_list, y, x, mode, V, fwhm, fmerit, step, inti, intf,

opt_npc = pclist[argm]
if verbose:
print 'Number of steps', len(pclist)
print('Number of steps', len(pclist))
msg = 'Optimal number of PCs = {}, for SNR={}'
print msg.format(opt_npc, snrlist[argm])
print
print(msg.format(opt_npc, snrlist[argm]))
print()
timing(start_time)

if full_output:
Expand All @@ -761,7 +761,7 @@ def grid(matrix, angle_list, y, x, mode, V, fwhm, fmerit, step, inti, intf,
ax2.set_ylabel('Flux in FWHM ap. [ADUs]')
ax2.minorticks_on()
ax2.grid('on', 'major', linestyle='solid', alpha=0.4)
print
print()

# automatic "clever" grid
else:
Expand Down Expand Up @@ -792,10 +792,9 @@ def grid(matrix, angle_list, y, x, mode, V, fwhm, fmerit, step, inti, intf,
ind = np.array(dfrsrd.index)

if verbose:
print 'Number of evaluated steps', ind.shape[0]
print('Number of evaluated steps', ind.shape[0])
msg = 'Optimal number of PCs = {}, for SNR={}'
print msg.format(opt_npc, snrlist3[argm])
print
print(msg.format(opt_npc, snrlist3[argm]), '\n')
timing(start_time)

if full_output:
Expand Down Expand Up @@ -832,7 +831,7 @@ def grid(matrix, angle_list, y, x, mode, V, fwhm, fmerit, step, inti, intf,
ax2.set_yscale('log')
ax2.grid('on', 'major', linestyle='solid', alpha=0.2)
#plt.savefig('figure.pdf', dpi=300, bbox_inches='tight')
print
print()

if mode == 'fullfr':
finalfr = pca(cube, angle_list, cube_ref=cube_ref, ncomp=opt_npc,
Expand Down Expand Up @@ -922,28 +921,28 @@ def pca_incremental(cubepath, angle_list=None, n=0, batch_size=None,

if verbose:
msg1 = "Cube with {:} frames ({:.3f} GB)"
print msg1.format(n_frames, hdulist[n].data.nbytes/1e9)
print(msg1.format(n_frames, hdulist[n].data.nbytes/1e9))
msg2 = "Batch size set to {:} frames ({:.3f} GB)"
print msg2.format(batch_size, hdulist[n].data[:batch_size].nbytes/1e9)
print
print(msg2.format(batch_size, hdulist[n].data[:batch_size].nbytes/1e9),
'\n')

res = n_frames%batch_size
for i in xrange(0, int(n_frames/batch_size)):
for i in range(0, int(n_frames/batch_size)):
intini = i*batch_size
intfin = (i+1)*batch_size
batch = hdulist[n].data[intini:intfin]
msg = 'Processing batch [{},{}] with shape {}'
if verbose:
print msg.format(intini, intfin, batch.shape)
print 'Batch size in memory = {:.3f} MB'.format(batch.nbytes/1e6)
print(msg.format(intini, intfin, batch.shape))
print('Batch size in memory = {:.3f} MB'.format(batch.nbytes/1e6))
matrix = prepare_matrix(batch, verbose=False)
ipca.partial_fit(matrix)
if res>0:
batch = hdulist[n].data[intfin:]
msg = 'Processing batch [{},{}] with shape {}'
if verbose:
print msg.format(intfin, n_frames, batch.shape)
print 'Batch size in memory = {:.3f} MB'.format(batch.nbytes/1e6)
print(msg.format(intfin, n_frames, batch.shape))
print('Batch size in memory = {:.3f} MB'.format(batch.nbytes/1e6))
matrix = prepare_matrix(batch, verbose=False)
ipca.partial_fit(matrix)

Expand All @@ -953,10 +952,9 @@ def pca_incremental(cubepath, angle_list=None, n=0, batch_size=None,
mean = ipca.mean_.reshape(batch.shape[1], batch.shape[2])

if verbose:
print
print 'Reconstructing and obtaining residuals'
print('\nReconstructing and obtaining residuals')
medians = []
for i in xrange(0, int(n_frames/batch_size)):
for i in range(0, int(n_frames/batch_size)):
intini = i*batch_size
intfin = (i+1)*batch_size
batch = hdulist[n].data[intini:intfin]
Expand Down
Loading

0 comments on commit b1db7c6

Please sign in to comment.