From de428fbfed57becbf83c0074f1a923545baf5348 Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 6 Mar 2019 12:27:10 -0500 Subject: [PATCH 01/33] trying out python 3.7 --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index b6c3194f0..dda75eb2b 100644 --- a/environment.yml +++ b/environment.yml @@ -1,7 +1,7 @@ channels: - conda-forge dependencies: -- python=3.6 +- python=3.7 - tensorflow # Must list before keras to preempt default keras backend - bokeh - cython From aa417e5977248393264fbab61a0fc66a8af25093 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Wed, 3 Apr 2019 12:19:04 -0400 Subject: [PATCH 02/33] type annotation effort 2019-04-03 --- caiman/base/movies.py | 179 +++++++++++++++++++------------- caiman/cluster.py | 52 ++++------ caiman/components_evaluation.py | 84 ++++----------- caiman/mmapping.py | 7 +- caiman/motion_correction.py | 73 +------------ 5 files changed, 149 insertions(+), 246 deletions(-) diff --git a/caiman/base/movies.py b/caiman/base/movies.py index 8eac3293f..98d004eaa 100644 --- a/caiman/base/movies.py +++ b/caiman/base/movies.py @@ -42,7 +42,7 @@ import sys import tifffile from tqdm import tqdm -from typing import List, Tuple +from typing import Any, Dict, List, Tuple, Union import warnings from zipfile import ZipFile @@ -226,9 +226,9 @@ def bin_median(self, window=10): num_frames = num_windows * window return np.nanmedian(np.nanmean(np.reshape(self[:num_frames], (window, num_windows, d1, d2)), axis=0), axis=0) - def extract_shifts(self, max_shift_w=5, max_shift_h=5, template=None, method='opencv'): + def extract_shifts(self, max_shift_w:int=5, max_shift_h:int=5, template=None, method:str='opencv') -> Tuple[List, List]: """ - Performs motion corretion using the opencv matchtemplate function. At every iteration a template is built by taking the median of all frames and then used to align the other frames. + Performs motion correction using the opencv matchtemplate function. At every iteration a template is built by taking the median of all frames and then used to align the other frames. Args: max_shift_w,max_shift_h: maximum pixel shifts allowed when correcting in the width and height direction @@ -442,7 +442,7 @@ def crop(self, crop_top=0, crop_bottom=0, crop_left=0, crop_right=0, crop_begin= t, h, w = self.shape return self[crop_begin:t - crop_end, crop_top:h - crop_bottom, crop_left:w - crop_right] - def computeDFF(self, secsWindow=5, quantilMin=8, method='only_baseline', order='F'): + def computeDFF(self, secsWindow:int=5, quantilMin:int=8, method:str='only_baseline', order:str='F') -> Tuple[Any, cm.movie]: """ compute the DFF of the movie or remove baseline @@ -513,7 +513,7 @@ def computeDFF(self, secsWindow=5, quantilMin=8, method='only_baseline', order=' logging.debug('Final Size Movie:' + np.str(self.shape)) return mov_out, movie(movBL, fr=self.fr, start_time=self.start_time, meta_data=self.meta_data, file_name=self.file_name) - def NonnegativeMatrixFactorization(self, n_components=30, init='nndsvd', beta=1, tol=5e-7, sparseness='components', **kwargs): + def NonnegativeMatrixFactorization(self, n_components:int=30, init:str='nndsvd', beta:int=1, tol=5e-7, sparseness:str='components', **kwargs) -> Tuple[np.ndarray, np.ndarray]: """ See documentation for scikit-learn NMF """ @@ -532,7 +532,7 @@ def NonnegativeMatrixFactorization(self, n_components=30, init='nndsvd', beta=1, return space_components, time_components - def online_NMF(self, n_components=30, method='nnsc', lambda1=100, iterations=-5, model=None, **kwargs): + def online_NMF(self, n_components:int=30, method:str='nnsc', lambda1:int=100, iterations:int=-5, model=None, **kwargs) -> Tuple[np.ndarray, np.ndarray]: """ Method performing online matrix factorization and using the spams (http://spams-devel.gforge.inria.fr/doc-python/html/index.html) package from Inria. @@ -585,7 +585,7 @@ def online_NMF(self, n_components=30, method='nnsc', lambda1=100, iterations=-5, return time_comps, np.array(space_comps) - def IPCA(self, components=50, batch=1000): + def IPCA(self, components:int=50, batch:int=1000) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ Iterative Principal Component analysis, see sklearn.decomposition.incremental_pca @@ -630,7 +630,7 @@ def IPCA(self, components=50, batch=1000): return eigenseries, eigenframes, proj_frame_vectors - def IPCA_stICA(self, componentsPCA=50, componentsICA=40, batch=1000, mu=1, ICAfun='logcosh', **kwargs): + def IPCA_stICA(self, componentsPCA:int=50, componentsICA:int=40, batch:int=1000, mu:float=1, ICAfun:str='logcosh', **kwargs) -> np.ndarray: """ Compute PCA + ICA a la Mukamel 2009. @@ -648,6 +648,7 @@ def IPCA_stICA(self, componentsPCA=50, componentsICA=40, batch=1000, mu=1, ICAfu Returns: ind_frames [components, height, width] = array of independent component "eigenframes" """ + # FIXME: mu defaults to 1, not 0.05 eigenseries, eigenframes, _proj = self.IPCA(componentsPCA, batch) # normalize the series @@ -675,7 +676,7 @@ def IPCA_stICA(self, componentsPCA=50, componentsICA=40, batch=1000, mu=1, ICAfu return ind_frames - def IPCA_denoise(self, components=50, batch=1000): + def IPCA_denoise(self, components:int=50, batch:int=1000): """ Create a denoised version of the movie using only the first 'components' components """ @@ -684,7 +685,7 @@ def IPCA_denoise(self, components=50, batch=1000): clean_vectors.T), np.shape(self)), **self.__dict__) return self - def IPCA_io(self, n_components=50, fun='logcosh', max_iter=1000, tol=1e-20): + def IPCA_io(self, n_components:int=50, fun:str='logcosh', max_iter:int=1000, tol=1e-20) -> np.ndarray: """ DO NOT USE STILL UNDER DEVELOPMENT """ pca_comp = n_components @@ -707,11 +708,11 @@ def IPCA_io(self, n_components=50, fun='logcosh', max_iter=1000, tol=1e-20): return mask - def local_correlations(self, eight_neighbours=False, swap_dim=True, frames_per_chunk=1500, order_mean=1): + def local_correlations(self, eight_neighbours:bool=False, swap_dim:bool=True, frames_per_chunk:int=1500, order_mean=1) -> np.ndarray: """Computes the correlation image for the input dataset Y Args: - Y: np.ndarray (3D or 4D) + self: np.ndarray (3D or 4D) Input movie data in 3D or 4D format eight_neighbours: Boolean @@ -722,6 +723,10 @@ def local_correlations(self, eight_neighbours=False, swap_dim=True, frames_per_c True indicates that time is listed in the last axis of Y (matlab format) and moves it in the front + frames_per_chunk: int (undocumented) + + order_mean: (undocumented) + Returns: rho: d1 x d2 [x d3] matrix, cross-correlation with adjacent pixels @@ -753,7 +758,7 @@ def local_correlations(self, eight_neighbours=False, swap_dim=True, frames_per_c return Cn - def partition_FOV_KMeans(self, tradeoff_weight=.5, fx=.25, fy=.25, n_clusters=4, max_iter=500): + def partition_FOV_KMeans(self, tradeoff_weight:float=.5, fx:float=.25, fy:float=.25, n_clusters:int=4, max_iter:int=500) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ Partition the FOV in clusters that are grouping pixels close in space and in mutual correlation @@ -767,7 +772,7 @@ def partition_FOV_KMeans(self, tradeoff_weight=.5, fx=.25, fy=.25, n_clusters=4, Returns: fovs:array 2D encoding the partitions of the FOV - mcoef: matric of pairwise correlation coefficients + mcoef: matrix of pairwise correlation coefficients distanceMatrix: matrix of picel distances """ @@ -790,7 +795,7 @@ def partition_FOV_KMeans(self, tradeoff_weight=.5, fx=.25, fy=.25, n_clusters=4, 1., fx), old_div(1., fy), interpolation=cv2.INTER_NEAREST) return np.uint8(fovs), mcoef, distanceMatrix - def extract_traces_from_masks(self, masks): + def extract_traces_from_masks(self, masks:np.ndarray) -> trace: """ Args: masks: array, 3D with each 2D slice bein a mask (integer or fractional) @@ -856,7 +861,7 @@ def resize(self, fx=1, fy=1, fz=1, interpolation=cv2.INTER_AREA): return self - def guided_filter_blur_2D(self, guide_filter, radius=5, eps=0): + def guided_filter_blur_2D(self, guide_filter, radius:int=5, eps=0): """ performs guided filtering on each frame. See opencv documentation of cv2.ximgproc.guidedFilter """ @@ -868,7 +873,7 @@ def guided_filter_blur_2D(self, guide_filter, radius=5, eps=0): return self - def bilateral_blur_2D(self, diameter=5, sigmaColor=10000, sigmaSpace=0): + def bilateral_blur_2D(self, diameter:int=5, sigmaColor:int=10000, sigmaSpace=0): """ performs bilateral filtering on each frame. See opencv documentation of cv2.bilateralFilter """ @@ -933,16 +938,12 @@ def median_blur_2D(self, kernel_size=3): return self - def resample(self): - # FIXME what is this function used for? - pass - - def to_2D(self, order='F'): + def to_2D(self, order='F') -> np.ndarray: [T, d1, d2] = self.shape d = d1 * d2 return np.reshape(self, (T, d), order=order) - def zproject(self, method='mean', cmap=pl.cm.gray, aspect='auto', **kwargs): + def zproject(self, method='mean', cmap=pl.cm.gray, aspect='auto', **kwargs) -> np.ndarray: """ Compute and plot projection across time: @@ -956,7 +957,8 @@ def zproject(self, method='mean', cmap=pl.cm.gray, aspect='auto', **kwargs): Raises: Exception 'Method not implemented' """ - # todo: todocument + # TODO: make the imshow optional + # TODO: todocument if method == 'mean': zp = np.mean(self, axis=0) elif method == 'median': @@ -968,18 +970,32 @@ def zproject(self, method='mean', cmap=pl.cm.gray, aspect='auto', **kwargs): pl.imshow(zp, cmap=cmap, aspect=aspect, **kwargs) return zp - def play(self, gain=1, fr=None, magnification=1, offset=0, interpolation=cv2.INTER_LINEAR, - backend='opencv', do_loop=False, bord_px=None, q_max=100, q_min = 0, plot_text = False): + def play(self, gain:float=1, fr=None, magnification=1, offset=0, interpolation=cv2.INTER_LINEAR, + backend:str='opencv', do_loop:bool=False, bord_px=None, q_max=100, q_min = 0, plot_text:bool=False) -> None: """ Play the movie using opencv Args: - gain: adjust movie brightness + gain: adjust movie brightness + + fr: framerate, playing speed if different from original (inter frame interval in seconds) + + magnification: (undocumented) + + offset: (undocumented) - frate : playing speed if different from original (inter frame interval in seconds) + interpolation: (undocumented) backend: 'pylab' or 'opencv', the latter much faster + do_loop: Whether to loop the video + + bord_px: (undocumented) + + q_max, q_min: (undocumented) + + plot_text: (undocumented) + Raises: Exception 'Unknown backend!' """ @@ -987,7 +1003,7 @@ def play(self, gain=1, fr=None, magnification=1, offset=0, interpolation=cv2.INT if backend == 'pylab': logging.warning('*** WARNING *** SPEED MIGHT BE LOW. USE opencv backend if available') - gain *= 1. + gain *= 1. # convert to float in case we were passed an int if q_max < 100: maxmov = np.nanpercentile(self[0:10], q_max) else: @@ -1089,11 +1105,11 @@ def animate(i): for i in range(10): cv2.waitKey(100) -def load(file_name, fr=30, start_time=0, meta_data=None, subindices=None, - shape=None, var_name_hdf5='mov', in_memory=False, is_behavior=False, - bottom=0, top=0, left=0, right=0, channel = None, outtype=np.float32): +def load(file_name, fr:float=30, start_time:float=0, meta_data:Dict=None, subindices=None, + shape:Tuple[int,int]=None, var_name_hdf5:str='mov', in_memory:bool=False, is_behavior:bool=False, + bottom=0, top=0, left=0, right=0, channel=None, outtype=np.float32) -> cm.movie: """ - load movie from file. SUpports a variety of formats. tif, hdf5, npy and memory mapped. Matlab is experimental. + load movie from file. Supports a variety of formats. tif, hdf5, npy and memory mapped. Matlab is experimental. Args: file_name: string @@ -1114,12 +1130,19 @@ def load(file_name, fr=30, start_time=0, meta_data=None, subindices=None, shape: tuple of two values dimension of the movie along x and y if loading from a two dimensional numpy array - num_frames_sub_idx: - when reading sbx format (experimental and unstable) - var_name_hdf5: str if loading from hdf5 name of the variable to load + in_memory: (undocumented) + + is_behavior: (undocumented) + + bottom,top,left,right: (undocumented) + + channel: (undocumented) + + outtype: The data type for the movie + Returns: mov: caiman.movie @@ -1196,7 +1219,7 @@ def load(file_name, fr=30, start_time=0, meta_data=None, subindices=None, height = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT)) cv_failed = False - dims = [length, height, width] + dims:Union[Tuple,List] = [length, height, width] if length == 0 or width == 0 or height == 0: #CV failed to load cv_failed = True if subindices is not None: @@ -1372,10 +1395,10 @@ def rgb2gray(rgb): return movie(input_arr.astype(outtype), fr=fr, start_time=start_time, file_name=os.path.split(file_name)[-1], meta_data=meta_data) -def load_movie_chain(file_list, fr=30, start_time=0, +def load_movie_chain(file_list:List[str], fr:float=30, start_time=0, meta_data=None, subindices=None, - bottom=0, top=0, left=0, right=0, z_top = 0, - z_bottom = 0, is3D = False, channel=None, outtype=np.float32): + bottom=0, top=0, left=0, right=0, z_top=0, + z_bottom=0, is3D:bool=False, channel=None, outtype=np.float32) -> cm.movie: """ load movies from list of file names Args: @@ -1420,47 +1443,47 @@ def load_movie_chain(file_list, fr=30, start_time=0, mov.append(m) return ts.concatenate(mov, axis=0) +#### +# TODO: Consider pulling these functions that work with .mat files into a separate file -def loadmat_sbx(filename): +def loadmat_sbx(filename:str): """ - this function should be called instead of direct spio.loadmat + this wrapper should be called instead of directly calling spio.loadmat - as it cures the problem of not properly recovering python dictionaries - from mat files. It calls the function check keys to cure all entries + It solves the problem of not properly recovering python dictionaries + from mat files. It calls the function check keys to fix all entries which are still mat-objects """ data_ = loadmat(filename, struct_as_record=False, squeeze_me=True) - return _check_keys(data_) - + _check_keys(data_) + return data_ -def _check_keys(dict): +def _check_keys(checkdict:Dict) -> None: """ - checks if entries in dictionary rare mat-objects. If yes todict is called to change them to nested dictionaries + checks if entries in dictionary rare mat-objects. If yes todict is called to change them to nested dictionaries. + Modifies its parameter in-place. """ - for key in dict: - if isinstance(dict[key], scipy.io.matlab.mio5_params.mat_struct): - dict[key] = _todict(dict[key]) - - return dict + for key in checkdict: + if isinstance(checkdict[key], scipy.io.matlab.mio5_params.mat_struct): + checkdict[key] = _todict(checkdict[key]) - -def _todict(matobj): +def _todict(matobj) -> Dict: """ A recursive function which constructs from matobjects nested dictionaries """ - dict = {} + ret = {} for strg in matobj._fieldnames: elem = matobj.__dict__[strg] if isinstance(elem, scipy.io.matlab.mio5_params.mat_struct): - dict[strg] = _todict(elem) + ret[strg] = _todict(elem) else: - dict[strg] = elem - return dict + ret[strg] = elem + return ret -def sbxread(filename, k=0, n_frames=np.inf): +def sbxread(filename:str, k:int=0, n_frames=np.inf) -> np.ndarray: """ Args: filename: str @@ -1508,11 +1531,13 @@ def sbxread(filename, k=0, n_frames=np.inf): return x.transpose([2, 1, 0]) -def sbxreadskip(filename, skip): +def sbxreadskip(filename:str, skip) -> np.ndarray: """ Args: filename: str filename should be full path excluding .sbx + + skip: (undocumented) """ # Check if contains .sbx and if so just truncate if '.sbx' in filename: @@ -1562,11 +1587,12 @@ def sbxreadskip(filename, skip): return x.transpose([2, 1, 0]) -def sbxshape(filename): +def sbxshape(filename:str) -> Tuple[int,int,int]: """ Args: filename should be full path excluding .sbx """ + # TODO: Document meaning of return values # Check if contains .sbx and if so just truncate if '.sbx' in filename: @@ -1594,20 +1620,21 @@ def sbxshape(filename): return x -def to_3D(mov2D, shape, order='F'): +def to_3D(mov2D:np.ndarray, shape:Tuple, order='F') -> np.ndarray: """ transform a vectorized movie into a 3D shape """ return np.reshape(mov2D, shape, order=order) -#%% -def from_zip_file_to_movie(zipfile_name, start_end = None): +def from_zip_file_to_movie(zipfile_name:str, start_end:Tuple=None) -> cm.movie: ''' + Read/convert a movie from a zipfile. - @param zipfile_name: - @param start_end: tuple - start and end frame to extract - @return: + Args: + zipfile_name: Name of zipfile to read + start_end: Start and end frame to extract + Returns: + cm.movie ''' mov:List = [] print('unzipping file into movie object') @@ -1617,7 +1644,7 @@ def from_zip_file_to_movie(zipfile_name, start_end = None): counter = 0 with ZipFile(zipfile_name) as archive: for idx, entry in enumerate(archive.infolist()): - if idx >= start_end[0] and idx < start_end[1]: + if idx >= start_end[0] and idx < start_end[1]: # type: ignore # FIXME: default value would not work with this with archive.open(entry) as file: if counter == 0: img = np.array(Image.open(file)) @@ -1634,13 +1661,17 @@ def from_zip_file_to_movie(zipfile_name, start_end = None): return cm.movie(mov[:counter]) -def from_zipfiles_to_movie_lists(zipfile_name, max_frames_per_movie=3000, binary = False): +def from_zipfiles_to_movie_lists(zipfile_name:str, max_frames_per_movie:int=3000, binary:bool=False) -> List[str]: ''' Transform zip file into set of tif movies - @param zipfile_name: - @param max_frames_per_movie: - @return: + Args: + zipfile_name: Name of zipfile to read + max_frames_per_movie: (undocumented) + binary: (undocumented) + Returns: + List of filenames in the movie list ''' + # TODO: Filter contents of zipfile for appropriate filenames, so we can have a README in there with ZipFile(zipfile_name) as archive: num_frames_total = len(archive.infolist()) diff --git a/caiman/cluster.py b/caiman/cluster.py index 6c6a22ff2..d83e56a15 100644 --- a/caiman/cluster.py +++ b/caiman/cluster.py @@ -34,27 +34,13 @@ import subprocess import sys import time +from typing import Any, Dict, List, Optional, Tuple, Union from .mmapping import load_memmap logger = logging.getLogger(__name__) -def get_patches_from_image(img, shapes, overlaps): - # todo todocument - d1, d2 = np.shape(img) - rf = np.divide(shapes, 2) - _, coords_2d = extract_patch_coordinates(d1, d2, rf=rf, stride=overlaps) - imgs = np.empty(coords_2d.shape[:2], dtype=np.object) - - for idx_0, count_0 in enumerate(coords_2d): - for idx_1, count_1 in enumerate(count_0): - imgs[idx_0, idx_1] = img[count_1[0], count_1[1]] - - return imgs, coords_2d -#%% - - -def extract_patch_coordinates(dims, rf, stride, border_pix=0, indeces=[slice(None)]*2): +def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[List[int],Tuple], border_pix:int=0, indeces=[slice(None)]*2) -> Tuple[List, List]: """ Partition the FOV in patches and return the indexed in 2D and 1D (flatten, order='F') formats @@ -119,30 +105,30 @@ def extract_patch_coordinates(dims, rf, stride, border_pix=0, indeces=[slice(Non for i, c in enumerate(coords_flat): assert len(c) == np.prod(shapes[i]) - return map(np.sort, coords_flat), shapes + return list(map(np.sort, coords_flat)), shapes #%% -def apply_to_patch(mmap_file, shape, dview, rf, stride, function, *args, **kwargs): +def apply_to_patch(mmap_file, shape:Tuple[Any,Any,Any], dview, rf, stride, function, *args, **kwargs) -> Tuple[List,Any,Tuple]: """ apply function to patches in parallel or not Args: - file_name: string - full path to an npy file (2D, pixels x time) containing the movie + mmap_file: Memory-mapped variable + Variable handle, first thing returned by load_memmap() shape: tuple of three elements dimensions of the original movie across y, x, and time + dview: ipyparallel view on client + if None + rf: int half-size of the square patch in pixel stride: int amount of overlap between patches - dview: ipyparallel view on client - if None - Returns: results @@ -150,6 +136,8 @@ def apply_to_patch(mmap_file, shape, dview, rf, stride, function, *args, **kwarg Exception 'Something went wrong' """ + # This function is (currently) only used in some demos and use_cases code. Unclear if it's still necessary, + # and it's been broken for awhile. (_, d1, d2) = shape if not np.isscalar(rf): @@ -165,7 +153,7 @@ def apply_to_patch(mmap_file, shape, dview, rf, stride, function, *args, **kwarg stride2 = stride idx_flat, idx_2d = extract_patch_coordinates( - d1, d2, rf=(rf1, rf2), stride=(stride1, stride2)) + (d1, d2), rf=(rf1, rf2), stride=(stride1, stride2)) shape_grid = tuple(np.ceil( (d1 * 1. / (rf1 * 2 - stride1), d2 * 1. / (rf2 * 2 - stride2))).astype(np.int)) @@ -179,9 +167,9 @@ def apply_to_patch(mmap_file, shape, dview, rf, stride, function, *args, **kwarg args_in = [] for id_f, id_2d in zip(idx_flat[:], idx_2d[:]): - args_in.append((mmap_file.filename, id_f, id_2d, function, args, kwargs)) + logger.debug("Flat index is of length " + str(len(idx_flat))) if dview is not None: try: @@ -223,7 +211,7 @@ def function_place_holder(args_in): #%% -def start_server(slurm_script=None, ipcluster="ipcluster", ncpus=None): +def start_server(slurm_script:str=None, ipcluster:str="ipcluster", ncpus:int=None) -> None: """ programmatically start the ipyparallel server @@ -268,7 +256,7 @@ def start_server(slurm_script=None, ipcluster="ipcluster", ncpus=None): c.close() sys.stdout.write("start_server: done\n") -def shell_source(script): +def shell_source(script:str) -> None: """ Run a source-style bash script, copy resulting env vars to current process. """ # XXX This function is weird and maybe not a good idea. People easily might expect # it to handle conditionals. Maybe just make them provide a key-value file @@ -289,14 +277,18 @@ def shell_source(script): os.environ.update(env) pipe.stdout.close() -def stop_server(ipcluster='ipcluster', pdir=None, profile=None, dview=None): +def stop_server(ipcluster:str='ipcluster', pdir:str=None, profile:str=None, dview=None) -> None: """ programmatically stops the ipyparallel server Args: ipcluster : str ipcluster binary file name; requires 4 path separators on Windows - Default: "ipcluster" + Default: "ipcluster"a + + pdir : Undocumented + profile: Undocumented + dview: Undocumented """ if 'multiprocessing' in str(type(dview)): @@ -359,7 +351,7 @@ def stop_server(ipcluster='ipcluster', pdir=None, profile=None, dview=None): #%% -def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=False): +def setup_cluster(backend:str='multiprocessing', n_processes:int=None, single_thread:bool=False) -> Tuple[Any, Any, Optional[int]]: """Setup and/or restart a parallel cluster. Args: backend: str diff --git a/caiman/components_evaluation.py b/caiman/components_evaluation.py index 73bd8fc23..24e7a75d7 100644 --- a/caiman/components_evaluation.py +++ b/caiman/components_evaluation.py @@ -18,7 +18,7 @@ import scipy from scipy.sparse import csc_matrix from scipy.stats import norm -from typing import Any, List +from typing import Any, List, Tuple, Union import warnings from caiman.paths import caiman_datadir @@ -29,44 +29,6 @@ except: pass -def estimate_noise_mode(traces, robust_std=False, use_mode_fast=False, return_all=False): - """ estimate the noise in the traces under assumption that signals are sparse and only positive. The last dimension should be time. - - """ - # todo todocument - if use_mode_fast: - md = mode_robust_fast(traces, axis=1) - else: - md = mode_robust(traces, axis=1) - - ff1 = traces - md[:, None] - # only consider values under the mode to determine the noise standard deviation - ff1 = -ff1 * (ff1 < 0) - if robust_std: - # compute 25 percentile - ff1 = np.sort(ff1, axis=1) - ff1[ff1 == 0] = np.nan - Ns = np.round(np.sum(ff1 > 0, 1) * .5) - iqr_h = np.zeros(traces.shape[0]) - - for idx, _ in enumerate(ff1): - iqr_h[idx] = ff1[idx, -Ns[idx]] - - # approximate standard deviation as iqr/1.349 - sd_r = 2 * iqr_h / 1.349 - - else: - Ns = np.sum(ff1 > 0, 1) - sd_r = np.sqrt(old_div(np.sum(ff1**2, 1), Ns)) - - if return_all: - return md, sd_r - else: - return sd_r -# - - -#%% try: profile except: @@ -74,7 +36,7 @@ def profile(a): return a @profile -def compute_event_exceptionality(traces, robust_std=False, N=5, use_mode_fast=False, sigma_factor=3.): +def compute_event_exceptionality(traces:np.ndarray, robust_std:bool=False, N:int=5, use_mode_fast:bool=False, sigma_factor:float=3.) -> Tuple[np.ndarray, np.ndarray, Any, Any]: """ Define a metric and order components according to the probability of some "exceptional events" (like a spike). @@ -162,7 +124,7 @@ def compute_event_exceptionality(traces, robust_std=False, N=5, use_mode_fast=Fa #%% -def find_activity_intervals(C, Npeaks=5, tB=-3, tA=10, thres=0.3): +def find_activity_intervals(C, Npeaks:int=5, tB=-3, tA=10, thres:float=0.3) -> List: # todo todocument import peakutils K, T = np.shape(C) @@ -192,7 +154,7 @@ def find_activity_intervals(C, Npeaks=5, tB=-3, tA=10, thres=0.3): #%% -def classify_components_ep(Y, A, C, b, f, Athresh=0.1, Npeaks=5, tB=-3, tA=10, thres=0.3): +def classify_components_ep(Y, A, C, b, f, Athresh=0.1, Npeaks=5, tB=-3, tA=10, thres=0.3) -> Tuple[np.ndarray, List]: # todo todocument K, _ = np.shape(C) @@ -238,7 +200,7 @@ def classify_components_ep(Y, A, C, b, f, Athresh=0.1, Npeaks=5, tB=-3, tA=10, t #%% -def evaluate_components_CNN(A, dims, gSig, model_name=os.path.join(caiman_datadir(), 'model', 'cnn_model'), patch_size=50, loaded_model=None, isGPU=False): +def evaluate_components_CNN(A, dims, gSig, model_name:str=os.path.join(caiman_datadir(), 'model', 'cnn_model'), patch_size:int=50, loaded_model=None, isGPU:bool=False) -> Tuple[Any, np.array]: """ evaluate component quality using a CNN network """ @@ -248,12 +210,8 @@ def evaluate_components_CNN(A, dims, gSig, model_name=os.path.join(caiman_datadi os.environ['CUDA_VISIBLE_DEVICES'] = '-1' -# try: os.environ["KERAS_BACKEND"] = "tensorflow" from keras.models import model_from_json -# except: -# logging.error('PROBLEM LOADING KERAS: cannot use classifier') - if loaded_model is None: if os.path.isfile(os.path.join(caiman_datadir(), model_name + ".json")): @@ -291,8 +249,8 @@ def evaluate_components_CNN(A, dims, gSig, model_name=os.path.join(caiman_datadi #%% -def evaluate_components(Y, traces, A, C, b, f, final_frate, remove_baseline=True, N=5, robust_std=False, - Athresh=0.1, Npeaks=5, thresh_C=0.3, sigma_factor=3.): +def evaluate_components(Y:np.ndarray, traces:np.ndarray, A, C, b, f, final_frate, remove_baseline:bool=True, N:int=5, robust_std:bool=False, + Athresh:float=0.1, Npeaks:int=5, thresh_C:float=0.3, sigma_factor:float=3.) -> Tuple[Any, Any, Any, Any, Any, Any]: """ Define a metric and order components according to the probability of some "exceptional events" (like a spike). Such probability is defined as the likeihood of observing the actual trace value over N samples given an estimated noise distribution. @@ -309,11 +267,13 @@ def evaluate_components(Y, traces, A, C, b, f, final_frate, remove_baseline=True Y: ndarray movie x,y,t + traces: ndarray + Fluorescence traces + A,C,b,f: various types outputs of cnmf - traces: ndarray - Fluorescence traces + final_frate: (undocumented) remove_baseline: bool whether to remove the baseline in a rolling fashion *(8 percentile) @@ -413,17 +373,11 @@ def evaluate_components(Y, traces, A, C, b, f, final_frate, remove_baseline=True return fitness_raw, fitness_delta, erfc_raw, erfc_delta, r_values, significant_samples - -#%% - -def grouper(n, iterable, fillvalue=None): +def grouper(n:int, iterable, fillvalue:bool=None): "grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx" args = [iter(iterable)] * n return itertools.zip_longest(*args, fillvalue=fillvalue) -#%% - - def evaluate_components_placeholder(params): import caiman as cm fname, traces, A, C, b, f, final_frate, remove_baseline, N, robust_std, Athresh, Npeaks, thresh_C = params @@ -435,12 +389,9 @@ def evaluate_components_placeholder(params): return fitness_raw, fitness_delta, [], [], r_values, significant_samples -#%% - - def estimate_components_quality_auto(Y, A, C, b, f, YrA, frate, decay_time, gSig, dims, dview=None, min_SNR=2, r_values_min=0.9, r_values_lowest=-1, Npeaks=10, use_cnn=True, thresh_cnn_min=0.95, thresh_cnn_lowest=0.1, - thresh_fitness_delta=-20., min_SNR_reject=0.5, gSig_range = None): + thresh_fitness_delta=-20., min_SNR_reject=0.5, gSig_range = None) -> Tuple[List, List, float, float, float]: ''' estimates the quality of component automatically Args: @@ -516,7 +467,7 @@ def estimate_components_quality_auto(Y, A, C, b, f, YrA, frate, decay_time, gSig traces = C + YrA - _, _, fitness_raw, _, r_values = estimate_components_quality( + _, _, fitness_raw, _, r_values = estimate_components_quality( # type: ignore # mypy cannot reason about return_all traces, Y, A, C, b, f, final_frate=frate, Npeaks=Npeaks, r_values_min=r_values_min, fitness_min=fitness_min, fitness_delta_min=thresh_fitness_delta, return_all=True, dview=dview, num_traces_per_group=50, N=N_samples) @@ -535,7 +486,7 @@ def select_components_from_metrics(A, dims, gSig, r_values, comp_SNR, min_SNR=2.5, min_SNR_reject=0.5, thresh_cnn_min=0.8, thresh_cnn_lowest=0.1, use_cnn=True, gSig_range=None, - neuron_class=1, predictions=None, **kwargs): + neuron_class=1, predictions=None, **kwargs) -> Tuple[Any, Any, Any]: '''Selects components based on pre-computed metrics. For each metric space correlation, trace SNR, and CNN classifier both an upper and a lower thresholds are considered. A component is accepted if and only if it @@ -585,8 +536,8 @@ def select_components_from_metrics(A, dims, gSig, r_values, comp_SNR, #%% def estimate_components_quality(traces, Y, A, C, b, f, final_frate=30, Npeaks=10, r_values_min=.95, - fitness_min=-100, fitness_delta_min=-100, return_all=False, N=5, - remove_baseline=True, dview=None, robust_std=False, Athresh=0.1, thresh_C=0.3, num_traces_per_group=20): + fitness_min=-100, fitness_delta_min=-100, return_all:bool=False, N=5, + remove_baseline=True, dview=None, robust_std=False, Athresh=0.1, thresh_C=0.3, num_traces_per_group=20) -> Union[Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.array], Tuple[np.ndarray, np.ndarray]]: """ Define a metric and order components according to the probability of some "exceptional events" (like a spike). Such probability is defined as the likeihood of observing the actual trace value over N samples given an estimated noise distribution. @@ -643,6 +594,7 @@ def estimate_components_quality(traces, Y, A, C, b, f, final_frate=30, Npeaks=10 float values representing correlation between component and spatial mask obtained by averaging important points """ + # TODO: Consider always returning it all and let the caller ignore what it does not want if 'memmap' not in str(type(Y)): logging.warning('NOT MEMORY MAPPED. FALLING BACK ON SINGLE CORE IMPLEMENTATION') diff --git a/caiman/mmapping.py b/caiman/mmapping.py index 51e7d982f..f9ee582c8 100644 --- a/caiman/mmapping.py +++ b/caiman/mmapping.py @@ -375,8 +375,7 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1 logging.debug('Distributing memory map over many files') # Here we make a bunch of memmap files in the right order. Same parameters - # TODO: Use separate variables to hold the list here vs the string we return - fname_new = cm.save_memmap_each(filenames, + fname_parts = cm.save_memmap_each(filenames, base_name = base_name, order = order, border_to_0 = border_to_0, @@ -388,14 +387,14 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1 slices = slices, add_to_movie = add_to_movie) else: - fname_new = filenames + fname_parts = filenames # The goal is to make a single large memmap file, which we do here if order == 'F': raise Exception('You cannot merge files in F order, they must be in C order for CaImAn') - fname_new = cm.save_memmap_join(fname_new, base_name=base_name, dview=dview, n_chunks=n_chunks) + fname_new = cm.save_memmap_join(fname_parts, base_name=base_name, dview=dview, n_chunks=n_chunks) else: # TODO: can be done online diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py index 967869905..bd6bb031c 100644 --- a/caiman/motion_correction.py +++ b/caiman/motion_correction.py @@ -40,7 +40,6 @@ """ from past.builtins import basestring -#%% from builtins import zip from builtins import map from builtins import str @@ -2474,76 +2473,6 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 extension = extension.lower() is_fiji = False -# if extension == '.tif' or extension == '.tiff': # check if tiff file -# with tifffile.TiffFile(fname) as tf: -# T = len(tf.pages) -# if T == 1: # Fiji-generated TIF -# is_fiji = True -# try: -# T, d1, d2 = tf[0].shape -# except: -# T, d1, d2 = tf.asarray().shape -# tf.close() -# else: -# d1, d2 = tf.pages[0].shape -# -# elif extension == '.sbx': # check if sbx file -# -# shape = cm.base.movies.sbxshape(name) -# d1 = shape[1] -# d2 = shape[0] -# T = shape[2] -# -# elif extension == '.sima': # check if sbx file -# import sima -# dataset = sima.ImagingDataset.load(fname) -# shape = dataset.sequences[0].shape -# d1 = shape[2] -# d2 = shape[3] -# T = shape[0] -# del dataset -# elif extension == '.npy': -# raise Exception('Numpy not supported at the moment') -# -# elif extension in ('.hdf5', '.h5'): -# with h5py.File(fname) as fl: -# fkeys = list(fl.keys()) -# if len(fkeys)==1: -# fsiz = fl[fkeys[0]].shape -# elif var_name_hdf5 in fkeys: -# fsiz = fl[var_name_hdf5].shape -# elif 'mov' in fkeys: -# fsiz = fl['mov'].shape -# elif 'imaging' in fkeys: -# fsiz = fl['imaging'].shape -# else: -# print(fkeys) -# raise Exception('Unsupported file key') -# if len(fsiz) == 3: -# T, d1, d2 = fsiz -# elif len(fsiz) == 5: -# T, _, d1, d2, _ = fsiz -# else: -# print(fsiz) -# raise Exception('Unsupported file shape') -# -# elif extension == '.avi': -# cap = cv2.VideoCapture(fname) -# try: -# T = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)) -# d2 = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) -# d1 = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)) -# except: -# logging.debug('Roll back top opencv 2') -# T = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_COUNT)) -# d2 = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_WIDTH)) -# d1 = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT)) -# cap.release() -# -# else: -# raise Exception( -# 'Unsupported file extension for parallel motion correction') - dims, T = cm.source_extraction.cnmf.utilities.get_file_size(fname, var_name_hdf5=var_name_hdf5) d1, d2 = dims @@ -2572,7 +2501,7 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 if save_movie: if base_name is None: base_name = os.path.split(fname)[1][:-4] - fname_tot = memmap_frames_filename(base_name, dims, T, order) + fname_tot:Optional[str] = memmap_frames_filename(base_name, dims, T, order) fname_tot = os.path.join(os.path.split(fname)[0], fname_tot) np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=prepare_shape(shape_mov), order=order) From b8c58d65c31c814134bea5292ef655fdb7c29e9b Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Wed, 3 Apr 2019 13:59:39 -0400 Subject: [PATCH 03/33] Correct spelling of indices across the codebase --- caiman/base/rois.py | 110 +++++++++--------- caiman/cluster.py | 8 +- caiman/source_extraction/cnmf/cnmf.py | 28 ++--- caiman/source_extraction/cnmf/estimates.py | 34 +++--- .../source_extraction/cnmf/initialization.py | 4 +- caiman/source_extraction/cnmf/map_reduce.py | 4 +- caiman/source_extraction/cnmf/merging.py | 10 +- caiman/source_extraction/cnmf/online_cnmf.py | 38 +++--- 8 files changed, 118 insertions(+), 118 deletions(-) diff --git a/caiman/base/rois.py b/caiman/base/rois.py index c1bb379f8..e14866c51 100644 --- a/caiman/base/rois.py +++ b/caiman/base/rois.py @@ -190,16 +190,16 @@ def nf_match_neurons_in_binary_masks(masks_gt, masks_comp, thresh_cost=.7, min_d Returns: idx_tp_1: - indeces true pos ground truth mask + indices true pos ground truth mask idx_tp_2: - indeces true pos comp + indices true pos comp idx_fn_1: - indeces false neg + indices false neg idx_fp_2: - indeces false pos + indices false pos """ @@ -357,16 +357,16 @@ def register_ROIs(A1, A2, dims, template1=None, template2=None, align_flag=True, Returns: matched_ROIs1: list - indeces of matched ROIs from session 1 + indices of matched ROIs from session 1 matched_ROIs2: list - indeces of matched ROIs from session 2 + indices of matched ROIs from session 2 non_matched1: list - indeces of non-matched ROIs from session 1 + indices of non-matched ROIs from session 1 non_matched2: list - indeces of non-matched ROIs from session 2 + indices of non-matched ROIs from session 2 performance: list (precision, recall, accuracy, f_1 score) with A1 taken as ground truth @@ -450,7 +450,7 @@ def register_ROIs(A1, A2, dims, template1=None, template2=None, align_flag=True, matches = matches[0] costs = costs[0] - #%% store indeces + #%% store indices idx_tp = np.where(np.array(costs) < thresh_cost)[0] if len(idx_tp) > 0: @@ -615,16 +615,16 @@ def register_multisession(A, dims, templates = [None], align_flag=True, return A_union, assignments, matchings -def extract_active_components(assignments, indeces, only = True): +def extract_active_components(assignments, indices, only = True): """ - Computes the indeces of components that were active in a specified set of + Computes the indices of components that were active in a specified set of sessions. Args: assignments: ndarray # of components X # of sessions assignments matrix returned by function register_multisession - indeces: list int + indices: list int set of sessions to look for active neurons. Session 1 corresponds to a pythonic index 0 etc @@ -635,14 +635,14 @@ def extract_active_components(assignments, indeces, only = True): Returns: components: list int - indeces of components + indices of components """ - components = np.where(np.isnan(assignments[:,indeces]).sum(-1) == 0)[0] + components = np.where(np.isnan(assignments[:,indices]).sum(-1) == 0)[0] if only: - not_inds = list(np.setdiff1d(range(assignments.shape[-1]), indeces)) + not_inds = list(np.setdiff1d(range(assignments.shape[-1]), indices)) not_comps = np.where(np.isnan(assignments[:,not_inds]).sum(-1) == len(not_inds))[0] components = np.intersect1d(components, not_comps) @@ -1240,10 +1240,10 @@ def detect_duplicates_and_subsets(binary_masks, predictions=None, r_values=None, logging.info(sz.shape) overlap = overlap/sz.T np.fill_diagonal(overlap, 0) - # pairs of duplicate indeces + # pairs of duplicate indices - indeces_orig = np.where((D < dist_thr) | ((overlap) >= thresh_subset)) - indeces_orig = [(a, b) for a, b in zip(indeces_orig[0], indeces_orig[1])] + indices_orig = np.where((D < dist_thr) | ((overlap) >= thresh_subset)) + indices_orig = [(a, b) for a, b in zip(indices_orig[0], indices_orig[1])] use_max_area = False if predictions is not None: @@ -1261,63 +1261,63 @@ def detect_duplicates_and_subsets(binary_masks, predictions=None, r_values=None, one, two = np.unravel_index(max_idx, overlap_tmp.shape) max_val = overlap_tmp[one, two] - indeces_to_keep:List = [] - indeces_to_remove = [] + indices_to_keep:List = [] + indices_to_remove = [] while max_val > 0: one, two = np.unravel_index(max_idx, overlap_tmp.shape) if metric[one] > metric[two]: - #indeces_to_keep.append(one) + #indices_to_keep.append(one) overlap_tmp[:,two] = 0 overlap_tmp[two,:] = 0 - indeces_to_remove.append(two) - #if two in indeces_to_keep: - # indeces_to_keep.remove(two) + indices_to_remove.append(two) + #if two in indices_to_keep: + # indices_to_keep.remove(two) else: overlap_tmp[:,one] = 0 overlap_tmp[one,:] = 0 - indeces_to_remove.append(one) - #indeces_to_keep.append(two) - #if one in indeces_to_keep: - # indeces_to_keep.remove(one) + indices_to_remove.append(one) + #indices_to_keep.append(two) + #if one in indices_to_keep: + # indices_to_keep.remove(one) max_idx = np.argmax(overlap_tmp) one, two = np.unravel_index(max_idx, overlap_tmp.shape) max_val = overlap_tmp[one, two] - #indeces_to_remove = np.setdiff1d(np.unique(indeces_orig),indeces_to_keep) - indeces_to_keep = np.setdiff1d(np.unique(indeces_orig),indeces_to_remove) + #indices_to_remove = np.setdiff1d(np.unique(indices_orig),indices_to_keep) + indices_to_keep = np.setdiff1d(np.unique(indices_orig),indices_to_remove) -# if len(indeces) > 0: +# if len(indices) > 0: # if use_max_area: # # if is to deal with tie breaks in case of same area -# indeces_keep = np.argmax([[overlap[sec, frst], overlap[frst, sec]] -# for frst, sec in indeces], 1) -# indeces_remove = np.argmin([[overlap[sec, frst], overlap[frst, sec]] -# for frst, sec in indeces], 1) +# indices_keep = np.argmax([[overlap[sec, frst], overlap[frst, sec]] +# for frst, sec in indices], 1) +# indices_remove = np.argmin([[overlap[sec, frst], overlap[frst, sec]] +# for frst, sec in indices], 1) # # # else: #use CNN -# indeces_keep = np.argmin([[metric[sec], metric[frst]] -# for frst, sec in indeces], 1) -# indeces_remove = np.argmax([[metric[sec], metric[frst]] -# for frst, sec in indeces], 1) +# indices_keep = np.argmin([[metric[sec], metric[frst]] +# for frst, sec in indices], 1) +# indices_remove = np.argmax([[metric[sec], metric[frst]] +# for frst, sec in indices], 1) # -# indeces_keep = np.unique([elms[ik] for ik, elms in -# zip(indeces_keep, indeces)]) -# indeces_remove = np.unique([elms[ik] for ik, elms in -# zip(indeces_remove, indeces)]) +# indices_keep = np.unique([elms[ik] for ik, elms in +# zip(indices_keep, indices)]) +# indices_remove = np.unique([elms[ik] for ik, elms in +# zip(indices_remove, indices)]) # -# multiple_appearance = np.intersect1d(indeces_keep,indeces_remove) +# multiple_appearance = np.intersect1d(indices_keep,indices_remove) # for mapp in multiple_appearance: -# indeces_remove.remove(mapp) +# indices_remove.remove(mapp) # else: -# indeces_keep = [] -# indeces_remove = [] -# indeces_keep = [] -# indeces_remove = [] +# indices_keep = [] +# indices_remove = [] +# indices_keep = [] +# indices_remove = [] - return indeces_orig, indeces_to_keep, indeces_to_remove, D, overlap + return indices_orig, indices_to_keep, indices_to_remove, D, overlap def detect_duplicates(file_name, dist_thr=0.1, FOV=(512, 512)): """ @@ -1331,9 +1331,9 @@ def detect_duplicates(file_name, dist_thr=0.1, FOV=(512, 512)): FOV: dimensions of the FOV Returns: - duplicates : list of indeces with duplicate entries + duplicates : list of indices with duplicate entries - ind_keep : list of kept indeces + ind_keep : list of kept indices """ rois = nf_read_roi_zip(file_name, FOV) @@ -1342,10 +1342,10 @@ def detect_duplicates(file_name, dist_thr=0.1, FOV=(512, 512)): np.reshape(rois, (rois.shape[0], np.prod(FOV))).T) D = distance_masks([sp_rois, sp_rois], [cm, cm], 10)[0] np.fill_diagonal(D, 1) - indeces = np.where(D < dist_thr) # pairs of duplicate indeces + indices = np.where(D < dist_thr) # pairs of duplicate indices - ind = list(np.unique(indeces[1][indeces[1] > indeces[0]])) + ind = list(np.unique(indices[1][indices[1] > indices[0]])) ind_keep = list(set(range(D.shape[0])) - set(ind)) - duplicates = list(np.unique(np.concatenate((indeces[0], indeces[1])))) + duplicates = list(np.unique(np.concatenate((indices[0], indices[1])))) return duplicates, ind_keep diff --git a/caiman/cluster.py b/caiman/cluster.py index d83e56a15..97b31f6b3 100644 --- a/caiman/cluster.py +++ b/caiman/cluster.py @@ -40,7 +40,7 @@ logger = logging.getLogger(__name__) -def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[List[int],Tuple], border_pix:int=0, indeces=[slice(None)]*2) -> Tuple[List, List]: +def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[List[int],Tuple], border_pix:int=0, indices=[slice(None)]*2) -> Tuple[List, List]: """ Partition the FOV in patches and return the indexed in 2D and 1D (flatten, order='F') formats @@ -56,9 +56,9 @@ def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[Lis degree of overlap of the patches """ - sl_start = [0 if sl.start is None else sl.start for sl in indeces] - sl_stop = [dim if sl.stop is None else sl.stop for (sl, dim) in zip(indeces, dims)] - sl_step = [1 for sl in indeces] # not used + sl_start = [0 if sl.start is None else sl.start for sl in indices] + sl_stop = [dim if sl.stop is None else sl.stop for (sl, dim) in zip(indices, dims)] + sl_step = [1 for sl in indices] # not used dims_large = dims dims = np.minimum(np.array(dims) - border_pix, sl_stop) - np.maximum(border_pix, sl_start) diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index e49676c12..f0bb70c11 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -295,7 +295,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p self.estimates = Estimates(A=Ain, C=Cin, b=b_in, f=f_in, dims=self.params.data['dims']) - def fit_file(self, motion_correct=False, indeces=[slice(None)]*2): + def fit_file(self, motion_correct=False, indices=[slice(None)]*2): fnames = self.params.get('data', 'fnames') if os.path.exists(fnames[0]): _, extension = os.path.splitext(fnames[0])[:2] @@ -331,7 +331,7 @@ def fit_file(self, motion_correct=False, indeces=[slice(None)]*2): images = np.reshape(Yr.T, [T] + list(dims), order='F') self.mmap_file = fname_new - return self.fit(images, indeces=indeces) + return self.fit(images, indices=indices) def refit(self, images, dview=None): """ @@ -353,7 +353,7 @@ def refit(self, images, dview=None): cnm.estimates = estimates return cnm.fit(images) - def fit(self, images, indeces=[slice(None), slice(None)]): + def fit(self, images, indices=[slice(None), slice(None)]): """ This method uses the cnmf algorithm to find sources in data. @@ -363,7 +363,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): Args: images : mapped np.ndarray of shape (t,x,y[,z]) containing the images that vary over time. - indeces: list of slice objects along dimensions (x,y[,z]) for processing only part of the FOV + indices: list of slice objects along dimensions (x,y[,z]) for processing only part of the FOV Returns: self: updated using the cnmf algorithm with C,A,S,b,f computed according to the given initial values @@ -378,16 +378,16 @@ def fit(self, images, indeces=[slice(None), slice(None)]): """ # Todo : to compartment - if isinstance(indeces, slice): - indeces = [indeces] - indeces = [slice(None)] + indeces - if len(indeces) < len(images.shape): - indeces = indeces + [slice(None)]*(len(images.shape) - len(indeces)) + if isinstance(indices, slice): + indices = [indices] + indices = [slice(None)] + indices + if len(indices) < len(images.shape): + indices = indices + [slice(None)]*(len(images.shape) - len(indices)) dims_orig = images.shape[1:] - dims_sliced = images[indeces].shape[1:] + dims_sliced = images[indices].shape[1:] is_sliced = (dims_orig != dims_sliced) if self.params.get('patch', 'rf') is None and (is_sliced or 'ndarray' in str(type(images))): - images = images[indeces] + images = images[indices] self.dview = None logging.warning("Parallel processing in a single patch " "is not available for loaded in memory or sliced" + @@ -507,7 +507,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): # embed in the whole FOV if is_sliced: FOV = np.zeros(dims_orig, order='C') - FOV[indeces[1:]] = 1 + FOV[indices[1:]] = 1 FOV = FOV.flatten(order='F') ind_nz = np.where(FOV>0)[0].tolist() self.estimates.A = self.estimates.A.tocsc() @@ -538,7 +538,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): gnb=self.params.get('init', 'nb'), border_pix=self.params.get('patch', 'border_pix'), low_rank_background=self.params.get('patch', 'low_rank_background'), del_duplicates=self.params.get('patch', 'del_duplicates'), - indeces=indeces) + indices=indices) self.estimates.bl, self.estimates.c1, self.estimates.g, self.estimates.neurons_sn = None, None, None, None logging.info("merging") @@ -608,7 +608,7 @@ def remove_components(self, ind_rm): Args: ind_rm : list - indeces of components to be removed + indices of components to be removed """ self.estimates.Ab, self.estimates.Ab_dense, self.estimates.CC, self.estimates.CY, self.M,\ diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index 949cfe152..2d221cb1d 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -80,10 +80,10 @@ def __init__(self, A=None, b=None, C=None, f=None, R=None, dims=None): contour plot for each spatial footprint idx_components: list - indeces of accepted components + indices of accepted components idx_components_bad: list - indeces of rejected components + indices of rejected components SNR_comp: np.ndarray trace SNR for each component @@ -669,10 +669,10 @@ def select_components(self, idx_components=None, use_object=False, save_discarde Args: idx_components: list - indeces of components to be kept + indices of components to be kept use_object: bool - Flag to use self.idx_components for reading the indeces. + Flag to use self.idx_components for reading the indices. save_discarded_components: bool whether to save the components from initialization so that they can be restored using the restore_discarded_components method @@ -779,7 +779,7 @@ class label for neuron shapes def evaluate_components(self, imgs, params, dview=None): """Computes the quality metrics for each component and stores the - indeces of the components that pass user specified thresholds. The + indices of the components that pass user specified thresholds. The various thresholds and parameters can be passed as inputs. If left empty then they are read from self.params.quality'] @@ -806,9 +806,9 @@ def evaluate_components(self, imgs, params, dview=None): Returns: self: esimates object self.idx_components: np.array - indeces of accepted components + indices of accepted components self.idx_components_bad: np.array - indeces of rejected components + indices of rejected components self.SNR_comp: np.array SNR values for each temporal trace self.r_values: np.array @@ -883,9 +883,9 @@ def filter_components(self, imgs, params, new_dict={}, dview=None): Returns: self: estimates object self.idx_components: np.array - indeces of accepted components + indices of accepted components self.idx_components_bad: np.array - indeces of rejected components + indices of rejected components self.SNR_comp: np.array SNR values for each temporal trace self.r_values: np.array @@ -918,9 +918,9 @@ def filter_components(self, imgs, params, new_dict={}, dview=None): return self def manual_merge(self, components, params): - ''' merge a given list of components. The indeces + ''' merge a given list of components. The indices of components are not pythonic, i.e., they start from 1. Moreover, - the indeces refer to the absolute indeces, i.e., the indeces before + the indices refer to the absolute indices, i.e., the indices before spliting the components in accepted and rejected. If you want to e.g. merge components 1 from idx_components and 10 from idx_components_bad you will to set @@ -1000,7 +1000,7 @@ def manual_merge(self, components, params): nr = self.C.shape[0] good_neurons = np.setdiff1d(list(range(nr)), neur_id) if self.idx_components is not None: - new_indeces = list(range(len(good_neurons), + new_indices = list(range(len(good_neurons), len(good_neurons) + nbmrg)) mapping_mat = np.zeros(nr) @@ -1010,7 +1010,7 @@ def manual_merge(self, components, params): new_idx_bad = [mapping_mat[i] for i in self.idx_components_bad if i in gn_] new_idx.sort() new_idx_bad.sort() - self.idx_components = np.array(new_idx + new_indeces, dtype=int) + self.idx_components = np.array(new_idx + new_indices, dtype=int) self.idx_components_bad = np.array(new_idx_bad, dtype=int) self.A = scipy.sparse.hstack((self.A.tocsc()[:, good_neurons], @@ -1090,7 +1090,7 @@ def remove_duplicates(self, predictions=None, r_values=None, dist_thr=0.1, min_d A_gt_thr_bin = (self.A_thr > 0).reshape([self.dims[0], self.dims[1], -1], order='F').transpose([2, 0, 1]) * 1. - duplicates_gt, indeces_keep_gt, indeces_remove_gt, D_gt, overlap_gt = detect_duplicates_and_subsets( + duplicates_gt, indices_keep_gt, indices_remove_gt, D_gt, overlap_gt = detect_duplicates_and_subsets( A_gt_thr_bin,predictions=predictions, r_values=r_values,dist_thr=dist_thr, min_dist=min_dist, thresh_subset=thresh_subset) print('Duplicates gt:' + str(len(duplicates_gt))) @@ -1101,14 +1101,14 @@ def remove_duplicates(self, predictions=None, r_values=None, dist_thr=0.1, min_d plt.imshow(A_gt_thr_bin[np.array(duplicates_gt).flatten()].sum(0)) plt.colorbar() plt.subplot(1, 3, 2) - plt.imshow(A_gt_thr_bin[np.array(indeces_keep_gt)[:]].sum(0)) + plt.imshow(A_gt_thr_bin[np.array(indices_keep_gt)[:]].sum(0)) plt.colorbar() plt.subplot(1, 3, 3) - plt.imshow(A_gt_thr_bin[np.array(indeces_remove_gt)[:]].sum(0)) + plt.imshow(A_gt_thr_bin[np.array(indices_remove_gt)[:]].sum(0)) plt.colorbar() plt.pause(1) - components_to_keep = np.delete(np.arange(self.A.shape[-1]), indeces_remove_gt) + components_to_keep = np.delete(np.arange(self.A.shape[-1]), indices_remove_gt) else: components_to_keep = np.arange(self.A.shape[-1]) diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index 265d6c307..1684d9709 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -644,9 +644,9 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, for s in ijSig], indexing='xy') arr = np.array([np.reshape(s, (1, np.size(s)), order='F').squeeze() for s in xySig], dtype=np.int) - indeces = np.ravel_multi_index(arr, d[0:-1], order='F') + indices = np.ravel_multi_index(arr, d[0:-1], order='F') - A[indeces, k] = np.reshape( + A[indices, k] = np.reshape( coef, (1, np.size(coef)), order='C').squeeze() Y[[slice(*a) for a in ijSig]] -= dataSig.copy() if k < nr - 1: diff --git a/caiman/source_extraction/cnmf/map_reduce.py b/caiman/source_extraction/cnmf/map_reduce.py index 719e24500..914c26a7b 100644 --- a/caiman/source_extraction/cnmf/map_reduce.py +++ b/caiman/source_extraction/cnmf/map_reduce.py @@ -144,7 +144,7 @@ def cnmf_patches(args_in): def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, border_pix=0, low_rank_background=True, - del_duplicates=False, indeces=[slice(None)]*3): + del_duplicates=False, indices=[slice(None)]*3): """Function that runs CNMF in patches Either in parallel or sequentially, and return the result for each. @@ -226,7 +226,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, params_copy.set('temporal', {'n_pixels_per_process': npx_per_proc}) idx_flat, idx_2d = extract_patch_coordinates( - dims, rfs, strides, border_pix=border_pix, indeces=indeces[1:]) + dims, rfs, strides, border_pix=border_pix, indices=indices[1:]) args_in = [] patch_centers = [] for id_f, id_2d in zip(idx_flat, idx_2d): diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index 5b162cd67..8da3d317e 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -134,12 +134,12 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, FF2 = A_corr > 0 C_corr = scipy.sparse.csc_matrix(A_corr.shape) for ii in range(nr): - overlap_indeces = A_corr[ii, :].nonzero()[1] - if len(overlap_indeces) > 0: + overlap_indices = A_corr[ii, :].nonzero()[1] + if len(overlap_indices) > 0: # we chesk the correlation of the calcium traces for eahc overlapping components corr_values = [scipy.stats.pearsonr(C[ii, :], C[jj, :])[ - 0] for jj in overlap_indeces] - C_corr[ii, overlap_indeces] = corr_values + 0] for jj in overlap_indices] + C_corr[ii, overlap_indices] = corr_values FF1 = (C_corr + C_corr.T) > thr FF3 = FF1.multiply(FF2) @@ -165,7 +165,7 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, # if not fast_merge: # Y_res = Y - A.dot(C) #residuals=background=noise if np.size(cor) > 1: - # we get the size (indeces) + # we get the size (indices) ind = np.argsort(np.squeeze(cor))[::-1] else: ind = [0] diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index a043a3785..b91a64514 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/caiman/source_extraction/cnmf/online_cnmf.py @@ -1466,34 +1466,34 @@ def get_candidate_components(sv, dims, Yres_buf, min_num_trial=3, gSig=(5, 5), ind = np.ravel_multi_index(ij, dims, order='C') ijSig = [[max(i - g, 0), min(i+g+1, d)] for i, g, d in zip(ij, gHalf, dims)] ijsig_all.append(ijSig) - indeces = np.ravel_multi_index(np.ix_(*[np.arange(ij[0], ij[1]) + indices = np.ravel_multi_index(np.ix_(*[np.arange(ij[0], ij[1]) for ij in ijSig]), dims, order='F').ravel(order='C') - # indeces_ = np.ravel_multi_index(np.ix_(*[np.arange(ij[0], ij[1]) + # indices_ = np.ravel_multi_index(np.ix_(*[np.arange(ij[0], ij[1]) # for ij in ijSig]), dims, order='C').ravel(order = 'C') - ain = np.maximum(mean_buff[indeces], 0) + ain = np.maximum(mean_buff[indices], 0) if sniper_mode: half_crop_cnn = tuple([int(np.minimum(gs*2, patch_size/2)) for gs in gSig]) ij_cnn = [min(max(ij_val,g_val),dim_val-g_val-1) for ij_val, g_val, dim_val in zip(ij,half_crop_cnn,dims)] ijSig_cnn = [[max(i - g, 0), min(i+g+1,d)] for i, g, d in zip(ij_cnn, half_crop_cnn, dims)] - indeces_cnn = np.ravel_multi_index(np.ix_(*[np.arange(ij[0], ij[1]) + indices_cnn = np.ravel_multi_index(np.ix_(*[np.arange(ij[0], ij[1]) for ij in ijSig_cnn]), dims, order='F').ravel(order = 'C') - ain_cnn = mean_buff[indeces_cnn] + ain_cnn = mean_buff[indices_cnn] else: compute_corr = True # determine when to compute corr coef na = ain.dot(ain) - # sv[indeces_] /= 1 # 0 + # sv[indices_] /= 1 # 0 if na: ain /= sqrt(na) Ain.append(ain) if compute_corr: - Y_patch.append(Yres_buf.T[indeces, :]) + Y_patch.append(Yres_buf.T[indices, :]) else: - all_indices.append(indeces) + all_indices.append(indices) idx.append(ind) if sniper_mode: Ain_cnn.append(ain_cnn) @@ -1589,7 +1589,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, ijSig = [[max(i - temp_g, 0), min(i + temp_g + 1, d)] for i, temp_g, d in zip(ij, gHalf, dims)] dims_ain = (np.abs(np.diff(ijSig[1])[0]), np.abs(np.diff(ijSig[0])[0])) - indeces = np.ravel_multi_index( + indices = np.ravel_multi_index( np.ix_(*[np.arange(ij[0], ij[1]) for ij in ijSig]), dims, order='F').ravel() @@ -1599,11 +1599,11 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, if Ab_dense is None: Ain = np.zeros((np.prod(dims), 1), dtype=np.float32) - Ain[indeces, :] = ain[:, None] + Ain[indices, :] = ain[:, None] ff = np.where((Ab.T.dot(Ain).T > thresh_overlap) [:, gnb:])[1] + gnb else: - ff = np.where(Ab_dense[indeces, gnb:M].T.dot( + ff = np.where(Ab_dense[indices, gnb:M].T.dot( ain).T > thresh_overlap)[0] + gnb if ff.size > 0: @@ -1620,7 +1620,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, # e.g. 1 * sigma * sqrt(1-sum(gamma)) corresponds roughly to the root mean square (non-zero) spike size, sqrt() # 2 * sigma * sqrt(1-sum(gamma)) corresponds roughly to the 95% percentile of (non-zero) spike sizes # 3 * sigma * sqrt(1-sum(gamma)) corresponds roughly to the 99.7% percentile of (non-zero) spike sizes - s_min = -s_min * sqrt((ain**2).dot(sn[indeces]**2)) * sqrt(1 - np.sum(g)) + s_min = -s_min * sqrt((ain**2).dot(sn[indices]**2)) * sqrt(1 - np.sum(g)) cin_res = cin_res.get_ordered() if accepted: @@ -1653,7 +1653,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, if not useOASIS: # lambda from Selesnick's 3*sigma*|K| rule # use noise estimate from init batch or use std_rr? - # sn_ = sqrt((ain**2).dot(sn[indeces]**2)) / sqrt(1 - g**2) + # sn_ = sqrt((ain**2).dot(sn[indices]**2)) / sqrt(1 - g**2) sn_ = std_rr oas = OASIS(np.ravel(g)[0], 3 * sn_ / (sqrt(1 - g**2) if np.size(g) == 1 else @@ -1667,12 +1667,12 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, oases.append(oas) - Ain_csc = csc_matrix((ain, (indeces, [0] * len(indeces))), (np.prod(dims), 1), dtype=np.float32) + Ain_csc = csc_matrix((ain, (indices, [0] * len(indices))), (np.prod(dims), 1), dtype=np.float32) if Ab_dense is None: groups = update_order(Ab, Ain, groups)[0] else: - groups = update_order(Ab_dense[indeces, :M], ain, groups)[0] - Ab_dense[indeces, M] = ain + groups = update_order(Ab_dense[indices, :M], ain, groups)[0] + Ab_dense[indices, M] = ain # faster version of scipy.sparse.hstack csc_append(Ab, Ain_csc) @@ -1684,7 +1684,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, Cf_ = Cf cin_circ_ = cin_circ - CY[M, indeces] = cin_.dot(Y_buf_[:, indeces]) / tt + CY[M, indices] = cin_.dot(Y_buf_[:, indices]) / tt # preallocate memory for speed up? CC1 = np.hstack([CC, Cf_.dot(cin_circ_ / tt)[:, None]]) @@ -1696,7 +1696,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, N = N + 1 M = M + 1 - Yres_buf[:, indeces] -= np.outer(cin, ain) + Yres_buf[:, indices] -= np.outer(cin, ain) # restrict blurring to region where component is located # update bigger region than neural patch to avoid boundary effects @@ -1765,7 +1765,7 @@ def remove_components_online(ind_rem, gnb, Ab, use_dense, Ab_dense, AtA, CY, Args: ind_rem list - indeces of components to be removed (starting from zero) + indices of components to be removed (starting from zero) gnb int number of global background components Ab csc_matrix From 869fd6ed8a6953edc21f9a78112e04461f525d81 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Wed, 3 Apr 2019 15:12:57 -0400 Subject: [PATCH 04/33] Adjust how to refer to movie class --- caiman/base/movies.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/caiman/base/movies.py b/caiman/base/movies.py index 98d004eaa..b9daad7bb 100644 --- a/caiman/base/movies.py +++ b/caiman/base/movies.py @@ -442,7 +442,7 @@ def crop(self, crop_top=0, crop_bottom=0, crop_left=0, crop_right=0, crop_begin= t, h, w = self.shape return self[crop_begin:t - crop_end, crop_top:h - crop_bottom, crop_left:w - crop_right] - def computeDFF(self, secsWindow:int=5, quantilMin:int=8, method:str='only_baseline', order:str='F') -> Tuple[Any, cm.movie]: + def computeDFF(self, secsWindow:int=5, quantilMin:int=8, method:str='only_baseline', order:str='F') -> Tuple[Any, movie]: """ compute the DFF of the movie or remove baseline @@ -1107,7 +1107,7 @@ def animate(i): def load(file_name, fr:float=30, start_time:float=0, meta_data:Dict=None, subindices=None, shape:Tuple[int,int]=None, var_name_hdf5:str='mov', in_memory:bool=False, is_behavior:bool=False, - bottom=0, top=0, left=0, right=0, channel=None, outtype=np.float32) -> cm.movie: + bottom=0, top=0, left=0, right=0, channel=None, outtype=np.float32) -> movie: """ load movie from file. Supports a variety of formats. tif, hdf5, npy and memory mapped. Matlab is experimental. @@ -1398,7 +1398,7 @@ def rgb2gray(rgb): def load_movie_chain(file_list:List[str], fr:float=30, start_time=0, meta_data=None, subindices=None, bottom=0, top=0, left=0, right=0, z_top=0, - z_bottom=0, is3D:bool=False, channel=None, outtype=np.float32) -> cm.movie: + z_bottom=0, is3D:bool=False, channel=None, outtype=np.float32) -> movie: """ load movies from list of file names Args: @@ -1414,7 +1414,7 @@ def load_movie_chain(file_list:List[str], fr:float=30, start_time=0, flag for 3d data (adds a fourth dimension) Returns: - movie: cm.movie + movie: movie movie corresponding to the concatenation og the input files """ @@ -1626,7 +1626,7 @@ def to_3D(mov2D:np.ndarray, shape:Tuple, order='F') -> np.ndarray: """ return np.reshape(mov2D, shape, order=order) -def from_zip_file_to_movie(zipfile_name:str, start_end:Tuple=None) -> cm.movie: +def from_zip_file_to_movie(zipfile_name:str, start_end:Tuple=None) -> movie: ''' Read/convert a movie from a zipfile. @@ -1634,7 +1634,7 @@ def from_zip_file_to_movie(zipfile_name:str, start_end:Tuple=None) -> cm.movie: zipfile_name: Name of zipfile to read start_end: Start and end frame to extract Returns: - cm.movie + movie ''' mov:List = [] print('unzipping file into movie object') From b44ce2872700d7a331091178ba6b1e7758df3a47 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Wed, 3 Apr 2019 15:25:17 -0400 Subject: [PATCH 05/33] Not sure how to fix type for base/movies so let's punt for now --- caiman/base/movies.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/caiman/base/movies.py b/caiman/base/movies.py index b9daad7bb..c901ba0fb 100644 --- a/caiman/base/movies.py +++ b/caiman/base/movies.py @@ -442,7 +442,7 @@ def crop(self, crop_top=0, crop_bottom=0, crop_left=0, crop_right=0, crop_begin= t, h, w = self.shape return self[crop_begin:t - crop_end, crop_top:h - crop_bottom, crop_left:w - crop_right] - def computeDFF(self, secsWindow:int=5, quantilMin:int=8, method:str='only_baseline', order:str='F') -> Tuple[Any, movie]: + def computeDFF(self, secsWindow:int=5, quantilMin:int=8, method:str='only_baseline', order:str='F') -> Tuple[Any, Any]: """ compute the DFF of the movie or remove baseline @@ -1107,7 +1107,7 @@ def animate(i): def load(file_name, fr:float=30, start_time:float=0, meta_data:Dict=None, subindices=None, shape:Tuple[int,int]=None, var_name_hdf5:str='mov', in_memory:bool=False, is_behavior:bool=False, - bottom=0, top=0, left=0, right=0, channel=None, outtype=np.float32) -> movie: + bottom=0, top=0, left=0, right=0, channel=None, outtype=np.float32) -> Any: """ load movie from file. Supports a variety of formats. tif, hdf5, npy and memory mapped. Matlab is experimental. @@ -1398,7 +1398,7 @@ def rgb2gray(rgb): def load_movie_chain(file_list:List[str], fr:float=30, start_time=0, meta_data=None, subindices=None, bottom=0, top=0, left=0, right=0, z_top=0, - z_bottom=0, is3D:bool=False, channel=None, outtype=np.float32) -> movie: + z_bottom=0, is3D:bool=False, channel=None, outtype=np.float32) -> Any: """ load movies from list of file names Args: @@ -1626,7 +1626,7 @@ def to_3D(mov2D:np.ndarray, shape:Tuple, order='F') -> np.ndarray: """ return np.reshape(mov2D, shape, order=order) -def from_zip_file_to_movie(zipfile_name:str, start_end:Tuple=None) -> movie: +def from_zip_file_to_movie(zipfile_name:str, start_end:Tuple=None) -> Any: ''' Read/convert a movie from a zipfile. From 52b34e52269219a5be965912c5965b7ebae83306 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Thu, 4 Apr 2019 15:50:35 -0400 Subject: [PATCH 06/33] More cleanups --- caiman/components_evaluation.py | 9 +++---- caiman/mmapping.py | 7 +++-- caiman/source_extraction/cnmf/estimates.py | 8 +++--- .../source_extraction/cnmf/initialization.py | 5 ++-- caiman/source_extraction/cnmf/map_reduce.py | 7 ++--- caiman/source_extraction/cnmf/merging.py | 6 ++--- caiman/source_extraction/cnmf/online_cnmf.py | 27 +++++++++---------- 7 files changed, 35 insertions(+), 34 deletions(-) diff --git a/caiman/components_evaluation.py b/caiman/components_evaluation.py index 24e7a75d7..b6a3cbab7 100644 --- a/caiman/components_evaluation.py +++ b/caiman/components_evaluation.py @@ -391,7 +391,7 @@ def evaluate_components_placeholder(params): def estimate_components_quality_auto(Y, A, C, b, f, YrA, frate, decay_time, gSig, dims, dview=None, min_SNR=2, r_values_min=0.9, r_values_lowest=-1, Npeaks=10, use_cnn=True, thresh_cnn_min=0.95, thresh_cnn_lowest=0.1, - thresh_fitness_delta=-20., min_SNR_reject=0.5, gSig_range = None) -> Tuple[List, List, float, float, float]: + thresh_fitness_delta=-20., min_SNR_reject=0.5, gSig_range = None) -> Tuple[np.array, np.array, float, float, float]: ''' estimates the quality of component automatically Args: @@ -480,13 +480,12 @@ def estimate_components_quality_auto(Y, A, C, b, f, YrA, frate, decay_time, gSig return idx_components, idx_components_bad, comp_SNR, r_values, cnn_values -#%% def select_components_from_metrics(A, dims, gSig, r_values, comp_SNR, r_values_min=0.8, r_values_lowest=-1, min_SNR=2.5, min_SNR_reject=0.5, thresh_cnn_min=0.8, thresh_cnn_lowest=0.1, use_cnn=True, gSig_range=None, - neuron_class=1, predictions=None, **kwargs) -> Tuple[Any, Any, Any]: + neuron_class=1, predictions=None, **kwargs) -> Tuple[np.array, np.array, Any]: '''Selects components based on pre-computed metrics. For each metric space correlation, trace SNR, and CNN classifier both an upper and a lower thresholds are considered. A component is accepted if and only if it @@ -533,11 +532,9 @@ def select_components_from_metrics(A, dims, gSig, r_values, comp_SNR, return idx_components.astype(np.int), idx_components_bad.astype(np.int), cnn_values -#%% - def estimate_components_quality(traces, Y, A, C, b, f, final_frate=30, Npeaks=10, r_values_min=.95, fitness_min=-100, fitness_delta_min=-100, return_all:bool=False, N=5, - remove_baseline=True, dview=None, robust_std=False, Athresh=0.1, thresh_C=0.3, num_traces_per_group=20) -> Union[Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.array], Tuple[np.ndarray, np.ndarray]]: + remove_baseline=True, dview=None, robust_std=False, Athresh=0.1, thresh_C=0.3, num_traces_per_group=20) -> Tuple[np.ndarray, ...]: """ Define a metric and order components according to the probability of some "exceptional events" (like a spike). Such probability is defined as the likeihood of observing the actual trace value over N samples given an estimated noise distribution. diff --git a/caiman/mmapping.py b/caiman/mmapping.py index f9ee582c8..2428b1dbd 100644 --- a/caiman/mmapping.py +++ b/caiman/mmapping.py @@ -34,7 +34,7 @@ def prepare_shape(mytuple:Tuple) -> Tuple: return tuple(map(lambda x: np.uint64(x), mytuple)) #%% -def load_memmap(filename:str, mode:str='r') -> Union[Tuple[Any, Tuple[int,int],int], Tuple[Any, Tuple[int,int,int], int]]: +def load_memmap(filename:str, mode:str='r') -> Tuple[Any,Tuple,int]: """ Load a memory mapped file created by the function save_memmap Args: @@ -68,7 +68,10 @@ def load_memmap(filename:str, mode:str='r') -> Union[Tuple[Any, Tuple[int,int],i ), int(fpart[-5]), int(fpart[-1]), fpart[-3] Yr = np.memmap(file_to_load, mode=mode, shape=prepare_shape(( d1 * d2 * d3, T)), dtype=np.float32, order=order) - return (Yr, (d1, d2), T) if d3 == 1 else (Yr, (d1, d2, d3), T) + if d3 == 1: + return (Yr, (d1, d2), T) + else: + return (Yr, (d1, d2, d3), T) else: logging.error("Unknown extension for file " + str(filename)) raise Exception('Unknown file extension (should be .mmap)') diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index 2d221cb1d..16d727746 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -6,11 +6,13 @@ @author: epnevmatikakis """ -import numpy as np +import logging import matplotlib.pyplot as plt +import numpy as np import scipy.sparse +from typing import List + import caiman -import logging from .utilities import detrend_df_f from .spatial import threshold_components from .merging import merge_iteration @@ -145,7 +147,7 @@ def __init__(self, A=None, b=None, C=None, f=None, R=None, dims=None): self.groups = None self.dims = dims - self.shifts = [] + self.shifts:List = [] self.A_thr = None self.discarded_components = None diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index 1684d9709..80eecf3de 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -33,6 +33,7 @@ from sklearn.decomposition import NMF, FastICA from sklearn.utils.extmath import randomized_svd, squared_norm import sys +from typing import List import caiman from caiman.source_extraction.cnmf.deconvolution import constrained_foopsi @@ -1608,8 +1609,8 @@ def get_indices_of_pixels_on_ring(pixel): downscale(b0.reshape(dims, order='F'), (ssub, ssub)).reshape((-1, 1), order='F') - indices = [] - data = [] + indices:List = [] + data:List = [] indptr = [0] for p in range(d1*d2): index = get_indices_of_pixels_on_ring(p) diff --git a/caiman/source_extraction/cnmf/map_reduce.py b/caiman/source_extraction/cnmf/map_reduce.py index 914c26a7b..d90fed568 100644 --- a/caiman/source_extraction/cnmf/map_reduce.py +++ b/caiman/source_extraction/cnmf/map_reduce.py @@ -22,12 +22,13 @@ from past.utils import old_div from copy import copy, deepcopy -from sklearn.decomposition import NMF import logging import numpy as np import os import scipy +from sklearn.decomposition import NMF import time +from typing import Set from ...mmapping import load_memmap from ...cluster import extract_patch_coordinates @@ -461,9 +462,9 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, # B_tot = scipy.sparse.coo_matrix(B_tot) F_tot *= nB[:, None] - processed_idx = set([]) + processed_idx:Set = set([]) # needed if a patch has more than 1 background component - processed_idx_prev = set([]) + processed_idx_prev:Set = set([]) for _b in np.arange(B_tot.shape[-1]): idx_mask = np.where(B_tot[:, _b])[0] idx_mask_repeat = processed_idx.intersection(idx_mask) diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index 8da3d317e..622726ab8 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -148,11 +148,11 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, FF3) # % extract connected components p = temporal_params['p'] - list_conxcomp = [] + list_conxcomp_initial = [] for i in range(nb): # we list them if np.sum(connected_comp == i) > 1: - list_conxcomp.append((connected_comp == i).T) - list_conxcomp = np.asarray(list_conxcomp).T + list_conxcomp_initial.append((connected_comp == i).T) + list_conxcomp = np.asarray(list_conxcomp_initial).T if list_conxcomp.ndim > 1: cor = np.zeros((np.shape(list_conxcomp)[1], 1)) diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index b91a64514..fdcdc0d19 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/caiman/source_extraction/cnmf/online_cnmf.py @@ -18,10 +18,9 @@ from builtins import range from builtins import str from builtins import zip -from math import sqrt -from time import time import cv2 +from math import sqrt import numpy as np from past.utils import old_div from scipy.ndimage import percentile_filter @@ -30,6 +29,8 @@ from scipy.stats import norm from sklearn.decomposition import NMF from sklearn.preprocessing import normalize +from time import time +from typing import List, Tuple import caiman # from caiman.source_extraction.cnmf import params @@ -203,12 +204,12 @@ def _prepare_object(self, Yr, T, new_dims=None, idx_components=None): min(self.params.get('online', 'init_batch'), self.params.get('online', 'minibatch_shape')) - 1), 0) self.estimates.groups = list(map(list, update_order(self.estimates.Ab)[0])) self.update_counter = 2**np.linspace(0, 1, self.N, dtype=np.float32) - self.time_neuron_added = [] + self.time_neuron_added:List = [] for nneeuu in range(self.N): self.time_neuron_added.append((nneeuu, self.params.get('online', 'init_batch'))) if self.params.get('online', 'dist_shape_update'): self.time_spend = 0 - self.comp_upd = [] + self.comp_upd:List = [] # setup per patch classifier if self.params.get('online', 'path_to_model') is None or self.params.get('online', 'sniper_mode') is False: @@ -694,12 +695,12 @@ def fit_online(self, **kwargs): extra_files = len(fls) - 1 init_files = 1 t = init_batch - self.Ab_epoch = [] + self.Ab_epoch:List = [] t_online = [] self.comp_upd = [] - self.t_shapes = [] - self.t_detect = [] - self.t_motion = [] + self.t_shapes:List = [] + self.t_detect:List = [] + self.t_motion:List = [] max_shifts_online = self.params.get('online', 'max_shifts_online') if extra_files == 0: # check whether there are any additional files process_files = fls[:init_files] # end processing at this file @@ -1416,8 +1417,8 @@ def get_candidate_components(sv, dims, Yres_buf, min_num_trial=3, gSig=(5, 5), idx = [] all_indices = [] ijsig_all = [] - cnn_pos = [] - local_maxima = [] + cnn_pos:List = [] + local_maxima:List = [] Y_patch = [] ksize = tuple([int(3 * i / 2) * 2 + 1 for i in gSig]) compute_corr = test_both @@ -1522,7 +1523,7 @@ def get_candidate_components(sv, dims, Yres_buf, min_num_trial=3, gSig=(5, 5), rval = corr(ain.copy(), np.mean(Ypx, -1)) if rval > rval_thr: keep_corr.append(i) - keep_final = list(set().union(keep_cnn, keep_corr)) + keep_final:List = list(set().union(keep_cnn, keep_corr)) if len(keep_final) > 0: Ain = np.stack(Ain)[keep_final] else: @@ -1813,15 +1814,11 @@ def initialize_movie_online(Y, K, gSig, rf, stride, base_name, Initialize movie using CNMF on minibatch. See CNMF parameters """ - _, d1, d2 = Y.shape - dims = (d1, d2) Yr = Y.to_2D().T # merging threshold, max correlation allowed # order of the autoregressive system - #T = Y.shape[0] base_name = base_name + '.mmap' fname_new = Y.save(base_name, order='C') - #% Yr, dims, T = caiman.load_memmap(fname_new) d1, d2 = dims images = np.reshape(Yr.T, [T] + list(dims), order='F') From 127bfbf026b5f73431cfbc4310c5dca7bc66e29a Mon Sep 17 00:00:00 2001 From: epnev Date: Thu, 11 Apr 2019 14:38:47 -0400 Subject: [PATCH 07/33] add some docstrings --- caiman/source_extraction/cnmf/cnmf.py | 37 ++++++++++++++++++++++++--- docs/source/Overview.rst | 4 +-- docs/source/core_functions.rst | 37 +++++++++++++++------------ 3 files changed, 56 insertions(+), 22 deletions(-) diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 9b1ac9e20..47aafd05d 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -298,6 +298,23 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p dims=self.params.data['dims']) def fit_file(self, motion_correct=False, indeces=[slice(None)]*2): + """ + This method packages the analysis pipeline (motion correction, memory + mapping, patch based CNMF processing) in a single method that can be + called on a specific (sequence of) file(s). It is assumed that the CNMF + object already contains a params object where the location of the files + and all the relevant parameters have been specified. The method does + not perform the quality evaluation step. Consult demo_pipeline for an + example. + + Args: + motion_correct (bool) + flag for performing motion correction + indeces (list of slice objects) + perform analysis only on a part of the FOV + Returns: + cnmf object with the current estimates + """ fnames = self.params.get('data', 'fnames') if os.path.exists(fnames[0]): _, extension = os.path.splitext(fnames[0])[:2] @@ -358,8 +375,7 @@ def refit(self, images, dview=None): def fit(self, images, indeces=[slice(None), slice(None)]): """ This method uses the cnmf algorithm to find sources in data. - - it is calling everyfunction from the cnmf folder + it is calling every function from the cnmf folder you can find out more at how the functions are called and how they are laid out at the ipython notebook Args: @@ -606,7 +622,8 @@ def save(self,filename): raise Exception("Filename not supported") def remove_components(self, ind_rm): - """remove a specified list of components from the OnACID CNMF object. + """ + Remove a specified list of components from the CNMF object. Args: ind_rm : list @@ -625,7 +642,10 @@ def remove_components(self, ind_rm): self.params.set('online', {'expected_comps': expected_comps}) def compute_residuals(self, Yr): - """compute residual for each component (variable YrA) + """ + Compute residual trace for each component (variable YrA). + WARNING: At the moment this method is valid only for the 2p processing + pipeline Args: Yr : np.ndarray @@ -909,6 +929,15 @@ def initialize(self, Y, **kwargs): return self def preprocess(self, Yr): + """ + Examines data to remove corrupted pixels and computes the noise level + estimate fo each pixel. + + Args: + Yr: np.array (or memmap.array) + 2d array of data (pixels x timesteps) typically in memory + mapped form + """ Yr, self.estimates.sn, self.estimates.g, self.estimates.psx = preprocess_data( Yr, dview=self.dview, **self.params.get_group('preprocess')) return Yr diff --git a/docs/source/Overview.rst b/docs/source/Overview.rst index ef789cadb..be0659199 100644 --- a/docs/source/Overview.rst +++ b/docs/source/Overview.rst @@ -2,7 +2,7 @@ Overview ========= .. image:: ../LOGOS/Caiman_logo_FI.png - :width: 200px + :width: 300px :align: right CaImAn is a Python toolbox for large scale **Ca**\ lcium **Im**\ aging data **An**\ alysis and behavioral analysis. @@ -34,8 +34,8 @@ Developers/Contributors CaImAn is being developed at the `Flatiron Institute `_ with numerous contributions from the broader community. The main developers are -* Andrea Giovannucci, University of North Carolina at Chapel Hill, previously at Flatiron Institute * Eftychios A. Pnevmatikakis, Flatiron Institute +* Andrea Giovannucci, University of North Carolina at Chapel Hill, previously at Flatiron Institute * Johannes Friedrich, Flatiron Institute * Pat Gunn, Flatiron Institute diff --git a/docs/source/core_functions.rst b/docs/source/core_functions.rst index fb574e822..e13e2b569 100644 --- a/docs/source/core_functions.rst +++ b/docs/source/core_functions.rst @@ -29,6 +29,8 @@ Functions that are required to operate the package at a basic level caiman.base.rois.register_multisession + caiman.source_extraction.cnmf.utilities.detrend_df_f + Movie Handling --------------- @@ -51,25 +53,15 @@ Timeseries Handling .. automethod:: timeseries.save .. autofunction:: concatenate -ROIs +ROIs --------------- .. currentmodule:: caiman.base.rois -.. automethod:: com -.. automethod:: extract_binary_masks_from_structural_channel -.. automethod:: register_ROIs -.. automethod:: register_multisession - - -Parallel Processing functions ------------------------------ - -.. currentmodule:: caiman.cluster - -.. autofunction:: apply_to_patch -.. autofunction:: start_server -.. autofunction:: stop_server +.. autofunction:: register_ROIs +.. autofunction:: register_multisession +.. autofunction:: com +.. autofunction:: extract_binary_masks_from_structural_channel Memory mapping @@ -160,6 +152,8 @@ CNMF .. automethod:: CNMF.deconvolve .. automethod:: CNMF.update_spatial .. automethod:: CNMF.update_temporal +.. automethod:: CNMF.compute_residuals +.. automethod:: CNMF.remove_components .. automethod:: CNMF.HALS4traces .. automethod:: CNMF.HALS4footprints .. automethod:: CNMF.merge_comps @@ -214,10 +208,21 @@ Merge components .. autofunction:: merge_components + Utilities --------------- .. currentmodule:: caiman.source_extraction.cnmf.utilities -.. autofunction:: detrend_df_f_auto +.. autofunction:: detrend_df_f .. autofunction:: update_order .. autofunction:: get_file_size + + +Parallel Processing functions +----------------------------- + +.. currentmodule:: caiman.cluster + +.. autofunction:: apply_to_patch +.. autofunction:: start_server +.. autofunction:: stop_server From 37bfb92025d9a02dd0da272ba931c6f36686bee7 Mon Sep 17 00:00:00 2001 From: epnev Date: Mon, 22 Apr 2019 13:05:18 -0400 Subject: [PATCH 08/33] tracking empty merged components --- caiman/source_extraction/cnmf/cnmf.py | 6 +++--- caiman/source_extraction/cnmf/merging.py | 3 ++- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 47aafd05d..bc31f0520 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -260,6 +260,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p # these are movie properties that will be refactored into the Movie object self.dims = None + self.empty_merged = None # these are member variables related to the CNMF workflow self.skip_refinement = skip_refinement @@ -585,7 +586,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): not_merged = np.setdiff1d(list(range(len(self.estimates.YrA))), np.unique(np.concatenate(self.estimates.merged_ROIs))) self.estimates.YrA = np.concatenate([self.estimates.YrA[not_merged], - np.array([self.estimates.YrA[m].mean(0) for m in self.estimates.merged_ROIs])]) + np.array([self.estimates.YrA[m].mean(0) for ind, m in enumerate(self.estimates.merged_ROIs) if not self.empty_merged[ind]])]) if self.params.get('init', 'nb') == 0: self.estimates.W, self.estimates.b0 = compute_W( Yr, self.estimates.A.toarray(), self.estimates.C, self.dims, @@ -597,7 +598,6 @@ def fit(self, images, indeces=[slice(None), slice(None)]): else: self.estimates.S = self.estimates.C else: - while len(self.estimates.merged_ROIs) > 0: self.merge_comps(Yr, mx=np.Inf) @@ -892,7 +892,7 @@ def merge_comps(self, Y, mx=50, fast_merge=True): """merges components """ self.estimates.A, self.estimates.C, self.estimates.nr, self.estimates.merged_ROIs, self.estimates.S, \ - self.estimates.bl, self.estimates.c1, self.estimates.neurons_sn, self.estimates.g=\ + self.estimates.bl, self.estimates.c1, self.estimates.neurons_sn, self.estimates.g, self.empty_merged=\ merge_components(Y, self.estimates.A, self.estimates.b, self.estimates.C, self.estimates.f, self.estimates.S, self.estimates.sn, self.params.get_group('temporal'), self.params.get_group('spatial'), dview=self.dview, diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index 5b162cd67..7ffd5e881 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -247,8 +247,9 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, else: logging.info('No more components merged!') merged_ROIs = [] + empty = [] - return A, C, nr, merged_ROIs, S, bl, c1, sn, g + return A, C, nr, merged_ROIs, S, bl, c1, sn, g, empty def merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx, indx, temporal_params): From 04758e56e7cfeda0dfc7e8ea28831eb7479c6c0d Mon Sep 17 00:00:00 2001 From: epnev Date: Mon, 22 Apr 2019 20:43:39 -0400 Subject: [PATCH 09/33] use sparse matrices during thresholding --- caiman/source_extraction/cnmf/spatial.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index b3473eb74..8ae56f88b 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -211,7 +211,8 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if i < np.prod(dims): pixel_groups.append([Y_name, C_name, sn, ind2_[i:np.prod(dims)], list( range(i, np.prod(dims))), method_ls, cct]) - A_ = np.zeros((d, nr + np.size(f, 0))) # init A_ + #A_ = np.zeros((d, nr + np.size(f, 0))) # init A_ + A_ = scipy.sparse.lil_matrix((d, nr + np.size(f, 0))) if dview is not None: if 'multiprocessing' in str(type(dview)): parallel_result = dview.map_async( @@ -230,7 +231,6 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, logging.info("thresholding components") A_ = threshold_components(A_, dims, dview=dview, medw=medw, thr_method=thr_method, maxthr=maxthr, nrgthr=nrgthr, extract_cc=extract_cc, se=se, ss=ss) - ff = np.where(np.sum(A_, axis=0) == 0) # remove empty components if np.size(ff) > 0: ff = ff[0] @@ -487,8 +487,7 @@ def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgth # dims and nm of neurones d, nr = np.shape(A) # instanciation of A thresh. - Ath = np.zeros((d, nr)) - + #Ath = np.zeros((d, nr)) pars = [] # fo each neurons A_1 = scipy.sparse.csc_matrix(A) @@ -505,11 +504,17 @@ def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgth else: res = list(map(threshold_components_parallel, pars)) - + res.sort(key = lambda x: x[1]) + indices = [] + indptr = [0] + data = [] for r in res: At, i = r - Ath[:, i] = At + indptr.append(indptr[-1]+At.indptr[-1]) + indices.extend(At.indices.tolist()) + data.extend(At.data.tolist()) + Ath = csc_matrix((data, indices, indptr), shape=(d, nr)) return Ath @@ -604,7 +609,7 @@ def threshold_components_parallel(pars): BW = BW.flatten() Ath2[BW] = Ath[BW] - return Ath2, i + return csr_matrix(Ath2), i # %% @@ -843,7 +848,7 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, d1, d2, d3 = dims d, nr = np.shape(A) A = csc_matrix(A) - dist_indicator = scipy.sparse.csc_matrix((d, nr),dtype= np.float32) + dist_indicator = scipy.sparse.lil_matrix((d, nr),dtype= np.float32) if method == 'ellipse': Coor = dict() @@ -927,7 +932,7 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, - return dist_indicator + return csc_matrix(dist_indicator) #%% def construct_dilate_parallel(pars): """ From 899b1d92a7ee53a4fff62cc8a0843909548f6bc4 Mon Sep 17 00:00:00 2001 From: epnev Date: Tue, 23 Apr 2019 09:51:43 -0400 Subject: [PATCH 10/33] fix indexing deprecation warning --- caiman/source_extraction/cnmf/cnmf.py | 4 ++-- caiman/source_extraction/cnmf/initialization.py | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 47aafd05d..5f63ced69 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -402,10 +402,10 @@ def fit(self, images, indeces=[slice(None), slice(None)]): if len(indeces) < len(images.shape): indeces = indeces + [slice(None)]*(len(images.shape) - len(indeces)) dims_orig = images.shape[1:] - dims_sliced = images[indeces].shape[1:] + dims_sliced = images[tuple(indeces)].shape[1:] is_sliced = (dims_orig != dims_sliced) if self.params.get('patch', 'rf') is None and (is_sliced or 'ndarray' in str(type(images))): - images = images[indeces] + images = images[tuple(indeces)] self.dview = None logging.warning("Parallel processing in a single patch " "is not available for loaded in memory or sliced" + diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index fe9b921c7..1587762bc 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -634,7 +634,7 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, for c in range(len(ij))] # we create an array of it (fl like) and compute the trace like the pixel ij trough time dataTemp = np.array( - Y[[slice(*a) for a in ijSig]].copy(), dtype=np.float32) + Y[tuple([slice(*a) for a in ijSig])].copy(), dtype=np.float32) traceTemp = np.array(np.squeeze(rho[ij]), dtype=np.float32) coef, score = finetune(dataTemp, traceTemp, nIter=nIter) @@ -649,25 +649,25 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, A[indeces, k] = np.reshape( coef, (1, np.size(coef)), order='C').squeeze() - Y[[slice(*a) for a in ijSig]] -= dataSig.copy() + Y[tuple([slice(*a) for a in ijSig])] -= dataSig.copy() if k < nr - 1: Mod = [[np.maximum(ij[c] - 2 * gHalf[c], 0), np.minimum(ij[c] + 2 * gHalf[c] + 1, d[c])] for c in range(len(ij))] ModLen = [m[1] - m[0] for m in Mod] Lag = [ijSig[c] - Mod[c][0] for c in range(len(ij))] dataTemp = np.zeros(ModLen) - dataTemp[[slice(*a) for a in Lag]] = coef + dataTemp[tuple([slice(*a) for a in Lag])] = coef dataTemp = imblur(dataTemp[..., np.newaxis], sig=gSig, siz=gSiz, kernel=kernel) temp = dataTemp * score.reshape([1] * (Y.ndim - 1) + [-1]) - rho[[slice(*a) for a in Mod]] -= temp.copy() + rho[tuple([slice(*a) for a in Mod])] -= temp.copy() if rolling_sum: rho_filt = scipy.signal.lfilter( - rolling_filter, 1., rho[[slice(*a) for a in Mod]]**2) - v[[slice(*a) for a in Mod]] = np.amax(rho_filt, axis=-1) + rolling_filter, 1., rho[tuple([slice(*a) for a in Mod])]**2) + v[tuple([slice(*a) for a in Mod])] = np.amax(rho_filt, axis=-1) else: - v[[slice(*a) for a in Mod]] = np.sum(rho[[slice(*a) - for a in Mod]]**2, axis=-1) + v[tuple([slice(*a) for a in Mod])] = \ + np.sum(rho[tuple([slice(*a) for a in Mod])]**2, axis=-1) res = np.reshape(Y, (np.prod(d[0:-1]), d[-1]), order='F') + med.flatten(order='F')[:, None] From 5f53fc7b53ccbdca93d84ab60b4e4f9a56403e64 Mon Sep 17 00:00:00 2001 From: epnev Date: Tue, 23 Apr 2019 10:45:11 -0400 Subject: [PATCH 11/33] checking for corrupted temporal components --- caiman/source_extraction/cnmf/temporal.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 86a7610e2..28e20479e 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -389,11 +389,17 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, if dview is not None and not('multiprocessing' in str(type(dview))): dview.results.clear() - if scipy.linalg.norm(Cin - C, 'fro') <= 1e-3*scipy.linalg.norm(C, 'fro'): - logging.info("stopping: overall temporal component not changing" + - " significantly") + try: + if scipy.linalg.norm(Cin - C, 'fro') <= 1e-3*scipy.linalg.norm(C, 'fro'): + logging.info("stopping: overall temporal component not changing" + + " significantly") + break + else: # we keep Cin and do the iteration once more + Cin = C.copy() + except ValueError: + logging.warning("Aborting updating of temporal components due" + + " to possible numerical issues.") + C = Cin.copy() break - else: # we keep Cin and do the iteration once more - Cin = C return C, S, bl, YrA, c1, sn, g, lam From 1db74dd187665816ea99d9f4e8a9274e5547ef70 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Tue, 23 Apr 2019 16:58:49 -0400 Subject: [PATCH 12/33] Recent versions of numpy changed the default value for allow_pickle in np.load() to False; this overrides the new default --- caiman/tests/comparison/comparison.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/caiman/tests/comparison/comparison.py b/caiman/tests/comparison/comparison.py index c8382df23..b95fed0ab 100644 --- a/caiman/tests/comparison/comparison.py +++ b/caiman/tests/comparison/comparison.py @@ -263,7 +263,7 @@ def save_with_compare(self, istruth=False, params=None, dview=None, Cn=None): else: # if not we create a comparison first try: - with np.load(file_path, encoding='latin1') as dt: + with np.load(file_path, encoding='latin1', allow_pickle=True) as dt: rig_shifts = dt['rig_shifts'][()] A_patch = dt['A_patch'][()] A_full = dt['A_full'][()] From ea700dba2fb0661895669d7d2f6b22494d3db80d Mon Sep 17 00:00:00 2001 From: epnev Date: Tue, 23 Apr 2019 21:11:31 -0400 Subject: [PATCH 13/33] faster construction of sparse A --- caiman/source_extraction/cnmf/spatial.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 8ae56f88b..062937237 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -208,11 +208,11 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, for i in range(0, np.prod(dims) - n_pixels_per_process + 1, n_pixels_per_process): pixel_groups.append([Y_name, C_name, sn, ind2_[i:i + n_pixels_per_process], list( range(i, i + n_pixels_per_process)), method_ls, cct, ]) - if i < np.prod(dims): - pixel_groups.append([Y_name, C_name, sn, ind2_[i:np.prod(dims)], list( - range(i, np.prod(dims))), method_ls, cct]) + if i + n_pixels_per_process < np.prod(dims): + pixel_groups.append([Y_name, C_name, sn, ind2_[(i + n_pixels_per_process):np.prod(dims)], list( + range(i + n_pixels_per_process, np.prod(dims))), method_ls, cct]) #A_ = np.zeros((d, nr + np.size(f, 0))) # init A_ - A_ = scipy.sparse.lil_matrix((d, nr + np.size(f, 0))) + #A_ = scipy.sparse.lil_matrix((d, nr + np.size(f, 0))) if dview is not None: if 'multiprocessing' in str(type(dview)): parallel_result = dview.map_async( @@ -223,10 +223,18 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, dview.results.clear() else: parallel_result = list(map(regression_ipyparallel, pixel_groups)) + data = [] + rows = [] + cols = [] for chunk in parallel_result: for pars in chunk: px, idxs_, a = pars - A_[px, idxs_] = a + #A_[px, idxs_] = a + nz = np.where(a>0)[0] + data.extend(a[nz]) + rows.extend(len(nz)*[px]) + cols.extend(idxs_[nz]) + A_ = scipy.sparse.coo_matrix((data, (rows, cols)), shape=(d, nr + np.size(f, 0))) logging.info("thresholding components") A_ = threshold_components(A_, dims, dview=dview, medw=medw, thr_method=thr_method, From f0b5f4c84238af75012a86f03e618291a3f0ba4c Mon Sep 17 00:00:00 2001 From: epnev Date: Tue, 23 Apr 2019 21:31:44 -0400 Subject: [PATCH 14/33] small optimizations --- caiman/source_extraction/cnmf/spatial.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 062937237..00daec74c 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -260,10 +260,10 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if 'memmap' in str(type(Y)): Y_resf = parallel_dot_product(Y, f.T, dview=dview, block_size=block_size_spat, num_blocks_per_run=num_blocks_per_run_spat) - \ - A_.dot(coo_matrix(C[:nr, :]).dot(f.T)) + A_.dot(C[:nr].dot(f.T)) else: # Y*f' - A*(C*f') - Y_resf = np.dot(Y, f.T) - A_.dot(coo_matrix(C[:nr, :]).dot(f.T)) + Y_resf = np.dot(Y, f.T) - A_.dot(C[:nr].dot(f.T)) if update_background_components: From 07ce88947103fec7e08296ae2aa9718a80a85cbd Mon Sep 17 00:00:00 2001 From: Rishi Bedi Date: Wed, 24 Apr 2019 08:41:45 +0000 Subject: [PATCH 15/33] NaN check in background component decomposition --- caiman/source_extraction/cnmf/map_reduce.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/caiman/source_extraction/cnmf/map_reduce.py b/caiman/source_extraction/cnmf/map_reduce.py index 719e24500..38e4125b7 100644 --- a/caiman/source_extraction/cnmf/map_reduce.py +++ b/caiman/source_extraction/cnmf/map_reduce.py @@ -427,7 +427,11 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, # np.random.rand(gnb - 1, T)] mdl = NMF(n_components=gnb, verbose=False, init='nndsvdar', tol=1e-10, max_iter=100, shuffle=False, random_state=1) + # Filter out nan components in the bg components + nan_components = np.any(np.isnan(F_tot), axis=1) + F_tot = F_tot[~nan_components, :] _ = mdl.fit_transform(F_tot).T + Bm = Bm[:, ~nan_components] f = mdl.components_.squeeze() f = np.atleast_2d(f) for _ in range(100): From d338564600e40785d9d3d4f3310ff7a5ae81b501 Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 24 Apr 2019 10:30:28 -0400 Subject: [PATCH 16/33] fixing bug when identifying empty thresholded components and tests --- caiman/source_extraction/cnmf/merging.py | 2 +- caiman/source_extraction/cnmf/spatial.py | 38 +++++++++++++++++------- caiman/tests/comparison/comparison.py | 4 +-- 3 files changed, 30 insertions(+), 14 deletions(-) diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index 5b162cd67..9a3fc6327 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -132,7 +132,7 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, A_corr.setdiag(0) A_corr = A_corr.tocsc() FF2 = A_corr > 0 - C_corr = scipy.sparse.csc_matrix(A_corr.shape) + C_corr = scipy.sparse.lil_matrix(A_corr.shape) for ii in range(nr): overlap_indeces = A_corr[ii, :].nonzero()[1] if len(overlap_indeces) > 0: diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 00daec74c..0c842d1c2 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -239,13 +239,14 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, logging.info("thresholding components") A_ = threshold_components(A_, dims, dview=dview, medw=medw, thr_method=thr_method, maxthr=maxthr, nrgthr=nrgthr, extract_cc=extract_cc, se=se, ss=ss) - ff = np.where(np.sum(A_, axis=0) == 0) # remove empty components + #ff = np.where(np.sum(A_, axis=0) == 0) # remove empty components + ff = np.asarray(A_.sum(0) == 0).nonzero()[1] if np.size(ff) > 0: - ff = ff[0] - logging.info('eliminating {0} empty spatial components'.format(len(ff))) - A_ = np.delete(A_, list(ff[ff < nr]), 1) - C = np.delete(C, list(ff[ff < nr]), 0) - nr = nr - len(ff[ff < nr]) + logging.info('removing {0} empty spatial component(s)'.format(ff.shape[0])) + if any(ff < nr): + A_ = np.delete(A_, list(ff[ff < nr]), 1) + C = np.delete(C, list(ff[ff < nr]), 0) + nr = nr - len(ff[ff < nr]) if update_background_components: if low_rank_background: background_ff = list(filter(lambda i: i >= nb, ff - nr)) @@ -254,6 +255,7 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, background_ff = list(filter(lambda i: i >= 0, ff - nr)) f = np.delete(f, background_ff, 0) b_in = np.delete(b_in, background_ff, 1) + A_ = A_[:, :nr] A_ = coo_matrix(A_) logging.info("Computing residuals") @@ -523,7 +525,6 @@ def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgth data.extend(At.data.tolist()) Ath = csc_matrix((data, indices, indptr), shape=(d, nr)) - return Ath @@ -599,7 +600,7 @@ def threshold_components_parallel(pars): # if we have deleted the element if BW.max() == 0: - return Ath2, i + return csr_matrix(Ath2), i # # we want to extract the largest connected component ( to remove small unconnected pixel ) if extract_cc: @@ -856,7 +857,8 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, d1, d2, d3 = dims d, nr = np.shape(A) A = csc_matrix(A) - dist_indicator = scipy.sparse.lil_matrix((d, nr),dtype= np.float32) +# dist_indicator = scipy.sparse.lil_matrix((d, nr),dtype= np.float32) +# dist_indicator = scipy.sparse.csc_matrix((d, nr), dtype=np.float32) if method == 'ellipse': Coor = dict() @@ -902,6 +904,9 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, dist_indicator = True * np.ones((d, nr)) elif method == 'dilate': + indptr = [0] + indices = [] + data = [] if dview is None: for i in range(nr): A_temp = np.reshape(A[:, i].toarray(), dims[::-1]) @@ -913,7 +918,12 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, else: A_temp = grey_dilation(A_temp, [1] * len(dims)) - dist_indicator[:, i] = scipy.sparse.coo_matrix(np.squeeze(np.reshape(A_temp, (d, 1)))[:,None] > 0) + nz = np.where(np.squeeze(np.reshape(A_temp, (d, 1)))[:, None] > 0)[0].tolist() + indptr.append(indptr[-1] + len(nz)) + indices.extend(nz) + data.extend(len(nz)*[True]) + #dist_indicator[:, i] = scipy.sparse.coo_matrix(np.squeeze(np.reshape(A_temp, (d, 1)))[:, None] > 0) + dist_indicator = csc_matrix((data, indices, indptr), shape=(d, nr)) else: logging.info('dilate in parallel...') @@ -931,8 +941,14 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, i = 0 for res in parallel_result: - dist_indicator[:, i] = res + #import pdb + #pdb.set_trace() + indptr.append(indptr[-1] + len(res.row)) + indices.extend(res.row) + data.extend(len(res.row)*[True]) + #dist_indicator[:, i] = res i += 1 + dist_indicator = csc_matrix((data, indices, indptr), shape=(d, nr)) else: raise Exception('Not implemented') diff --git a/caiman/tests/comparison/comparison.py b/caiman/tests/comparison/comparison.py index b95fed0ab..5552aab14 100644 --- a/caiman/tests/comparison/comparison.py +++ b/caiman/tests/comparison/comparison.py @@ -444,8 +444,8 @@ def cnmf(Cn, A_gt, A_test, C_gt, C_test, dims_gt, dims_test, dview=None, sensiti se=None, ss=None, dview=dview) # compute C using this A thr - A_test_thr = A_test_thr > 0 - A_gt_thr = A_gt_thr > 0 + A_test_thr = A_test_thr.toarray() > 0 + A_gt_thr = A_gt_thr.toarray() > 0 # we do not compute a threshold on the size of neurons C_test_thr = C_test C_gt_thr = C_gt From 1e397281d773a52114ef280995991f0dcf629edd Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Wed, 24 Apr 2019 11:00:57 -0400 Subject: [PATCH 17/33] Jenkins/Win: Handle any version of the visual studio linkage --- Jenkinsfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Jenkinsfile b/Jenkinsfile index 8524e011b..3ceb6af71 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -72,7 +72,7 @@ pipeline { steps { bat '%ANACONDA3%\\scripts\\conda info' bat '%ANACONDA3%\\scripts\\conda env create -q -f environment.yml -p %CONDA_ENV%' - bat 'del "%CONDA_ENV%\\etc\\conda\\activate.d\\vs2015_compiler_vars.bat' + bat 'del "%CONDA_ENV%\\etc\\conda\\activate.d\\vs*_compiler_vars.bat"' bat '%ANACONDA3%\\scripts\\activate %CONDA_ENV% && set KERAS_BACKEND=tensorflow && pip install . && copy caimanmanager.py %TEMP% && cd %TEMP% && set "CAIMAN_DATA=%TEMP%\\caiman_data" && (if exist caiman_data (rmdir caiman_data /s /q && echo "Removed old caiman_data" ) else (echo "Host is fresh")) && python caimanmanager.py install --force && python caimanmanager.py test' } } From 71156b0c8142e095d6b429897997e0106fd2dfa6 Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 24 Apr 2019 12:59:04 -0400 Subject: [PATCH 18/33] add method for removing cols from csc_matrix --- caiman/source_extraction/cnmf/spatial.py | 6 ++--- caiman/utils/utils.py | 34 ++++++++++++++++++++++-- 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 0c842d1c2..55849b1fc 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -34,6 +34,7 @@ from typing import List from ...mmapping import load_memmap, parallel_dot_product +from ...utils.utils import csc_column_remove def basis_denoising(y, c, boh, sn, id2_, px): @@ -244,7 +245,8 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if np.size(ff) > 0: logging.info('removing {0} empty spatial component(s)'.format(ff.shape[0])) if any(ff < nr): - A_ = np.delete(A_, list(ff[ff < nr]), 1) + #A_ = np.delete(A_, list(ff[ff < nr]), 1) + A_ = csc_column_remove(A_, list(ff[ff < nr])) C = np.delete(C, list(ff[ff < nr]), 0) nr = nr - len(ff[ff < nr]) if update_background_components: @@ -941,8 +943,6 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, i = 0 for res in parallel_result: - #import pdb - #pdb.set_trace() indptr.append(indptr[-1] + len(res.row)) indices.extend(res.row) data.extend(len(res.row)*[True]) diff --git a/caiman/utils/utils.py b/caiman/utils/utils.py index b7eb41ee9..4ea2f569b 100644 --- a/caiman/utils/utils.py +++ b/caiman/utils/utils.py @@ -30,7 +30,7 @@ import scipy from scipy.ndimage.filters import gaussian_filter from tifffile import TiffFile -from typing import Any, Dict, List, Tuple, Union +from typing import Any, Dict, List, Tuple, Union, Iterable try: cv2.setNumThreads(0) @@ -493,4 +493,34 @@ def recursively_load_dict_contents_from_group(h5file:h5py.File, path:str) -> Dic indptr[:]), shape[:]) else: ans[key] = recursively_load_dict_contents_from_group(h5file, path + key + '/') - return ans \ No newline at end of file + return ans + + +def csc_column_remove(A: scipy.sparse.csc_matrix, ind: Iterable[int]) -> scipy.sparse.csc_matrix: + """ Removes specified columns for a scipy.sparse csc_matrix + Args: + A: scipy.sparse.csc_matrix + Input matrix + ind: iterable[int] + list or np.array with columns to be removed + """ + d1, d2 = A.shape + if 'csc_matrix' not in str(type(A)): + logging.warning("Original matrix not in csc_format. Converting it" + + " anyway.") + A = scipy.sparse.csc_matrix(A) + indptr = A.indptr + ind_diff = np.diff(A.indptr).tolist() + ind_sort = sorted(ind, reverse=True) + data_list = [A.data[indptr[i]:indptr[i+1]] for i in range(d2)] + indices_list = [A.indices[indptr[i]:indptr[i+1]] for i in range(d2)] + for i in ind_sort: + del data_list[i] + del indices_list[i] + del ind_diff[i] + indptr_final = np.cumsum([0] + ind_diff) + data_final = [item for sublist in data_list for item in sublist] + indices_final = [item for sublist in indices_list for item in sublist] + A = scipy.sparse.csc_matrix((data_final, indices_final, indptr_final), + shape=[d1, d2 - len(ind)]) + return A From 4d631181a0e5926f9798d4dc4e516695f6b59ec3 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Wed, 24 Apr 2019 13:41:10 -0400 Subject: [PATCH 19/33] Jenkins/Windows: Only delete VisualStudio glue if it exists (py37 compat) --- Jenkinsfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Jenkinsfile b/Jenkinsfile index 3ceb6af71..4fe116dcc 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -72,7 +72,7 @@ pipeline { steps { bat '%ANACONDA3%\\scripts\\conda info' bat '%ANACONDA3%\\scripts\\conda env create -q -f environment.yml -p %CONDA_ENV%' - bat 'del "%CONDA_ENV%\\etc\\conda\\activate.d\\vs*_compiler_vars.bat"' + bat 'if exist "%CONDA_ENV%\\etc\\conda\\activate.d\\vs*_compiler_vars.bat" del "%CONDA_ENV%\\etc\\conda\\activate.d\\vs*_compiler_vars.bat"' bat '%ANACONDA3%\\scripts\\activate %CONDA_ENV% && set KERAS_BACKEND=tensorflow && pip install . && copy caimanmanager.py %TEMP% && cd %TEMP% && set "CAIMAN_DATA=%TEMP%\\caiman_data" && (if exist caiman_data (rmdir caiman_data /s /q && echo "Removed old caiman_data" ) else (echo "Host is fresh")) && python caimanmanager.py install --force && python caimanmanager.py test' } } From d5caf4288f6bcf9aa96a2ffeb0d26ecd4c87709c Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 24 Apr 2019 14:02:25 -0400 Subject: [PATCH 20/33] fixing a cyclical import --- caiman/source_extraction/cnmf/spatial.py | 4 +- caiman/utils/stats.py | 61 +++++++++++++++++------- caiman/utils/utils.py | 39 ++------------- 3 files changed, 51 insertions(+), 53 deletions(-) diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 55849b1fc..9547c6cb4 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -34,7 +34,7 @@ from typing import List from ...mmapping import load_memmap, parallel_dot_product -from ...utils.utils import csc_column_remove +from ...utils.stats import csc_column_remove def basis_denoising(y, c, boh, sn, id2_, px): @@ -516,7 +516,7 @@ def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgth else: res = list(map(threshold_components_parallel, pars)) - res.sort(key = lambda x: x[1]) + res.sort(key=lambda x: x[1]) indices = [] indptr = [0] data = [] diff --git a/caiman/utils/stats.py b/caiman/utils/stats.py index 9f0845ee7..6dd633932 100644 --- a/caiman/utils/stats.py +++ b/caiman/utils/stats.py @@ -15,6 +15,7 @@ pass import numpy as np +import scipy #%% @@ -182,15 +183,10 @@ def fnc(x): return df_percentile(x) Updated 1-23-2013 """ - -import scipy as sci -import scipy.optimize -import scipy.fftpack - def kde(data, N=None, MIN=None, MAX=None): # Parameters to set up the mesh on which to calculate - N = 2**12 if N is None else int(2**sci.ceil(sci.log2(N))) + N = 2**12 if N is None else int(2**scipy.ceil(scipy.log2(N))) if MIN is None or MAX is None: minimum = min(data) maximum = max(data) @@ -203,7 +199,7 @@ def kde(data, N=None, MIN=None, MAX=None): # Histogram the data to get a crude first approximation of the density M = len(data) - DataHist, bins = sci.histogram(data, bins=N, range=(MIN, MAX)) + DataHist, bins = scipy.histogram(data, bins=N, range=(MIN, MAX)) DataHist = DataHist/M DCTData = scipy.fftpack.dct(DataHist, norm=None) @@ -220,26 +216,57 @@ def kde(data, N=None, MIN=None, MAX=None): return None # Smooth the DCTransformed data using t_star - SmDCTData = DCTData*sci.exp(-sci.arange(N)**2*sci.pi**2*t_star/2) + SmDCTData = DCTData*scipy.exp(-scipy.arange(N)**2*scipy.pi**2*t_star/2) # Inverse DCT to get density density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R mesh = [(bins[i]+bins[i+1])/2 for i in range(N)] - bandwidth = sci.sqrt(t_star)*R + bandwidth = scipy.sqrt(t_star)*R - density = density/sci.trapz(density, mesh) + density = density/scipy.trapz(density, mesh) cdf = np.cumsum(density)*(mesh[1]-mesh[0]) return bandwidth, mesh, density, cdf + def fixed_point(t, M, I, a2): l=7 - I = sci.float64(I) - M = sci.float64(M) - a2 = sci.float64(a2) - f = 2*sci.pi**(2*l)*sci.sum(I**l*a2*sci.exp(-I*sci.pi**2*t)) + I = scipy.float64(I) + M = scipy.float64(M) + a2 = scipy.float64(a2) + f = 2*scipy.pi**(2*l)*scipy.sum(I**l*a2*scipy.exp(-I*scipy.pi**2*t)) for s in range(l, 1, -1): - K0 = sci.prod(range(1, 2*s, 2))/sci.sqrt(2*sci.pi) + K0 = scipy.prod(range(1, 2*s, 2))/scipy.sqrt(2*scipy.pi) const = (1 + (1/2)**(s + 1/2))/3 time=(2*const*K0/M/f)**(2/(3+2*s)) - f=2*sci.pi**(2*s)*sci.sum(I**s*a2*sci.exp(-I*sci.pi**2*time)) - return t-(2*M*sci.sqrt(sci.pi)*f)**(-2/5) \ No newline at end of file + f=2*scipy.pi**(2*s)*scipy.sum(I**s*a2*scipy.exp(-I*scipy.pi**2*time)) + return t-(2*M*scipy.sqrt(scipy.pi)*f)**(-2/5) + + +def csc_column_remove(A, ind): + """ Removes specified columns for a scipy.sparse csc_matrix + Args: + A: scipy.sparse.csc_matrix + Input matrix + ind: iterable[int] + list or np.array with columns to be removed + """ + d1, d2 = A.shape + if 'csc_matrix' not in str(type(A)): + logging.warning("Original matrix not in csc_format. Converting it" + + " anyway.") + A = scipy.sparse.csc_matrix(A) + indptr = A.indptr + ind_diff = np.diff(A.indptr).tolist() + ind_sort = sorted(ind, reverse=True) + data_list = [A.data[indptr[i]:indptr[i+1]] for i in range(d2)] + indices_list = [A.indices[indptr[i]:indptr[i+1]] for i in range(d2)] + for i in ind_sort: + del data_list[i] + del indices_list[i] + del ind_diff[i] + indptr_final = np.cumsum([0] + ind_diff) + data_final = [item for sublist in data_list for item in sublist] + indices_final = [item for sublist in indices_list for item in sublist] + A = scipy.sparse.csc_matrix((data_final, indices_final, indptr_final), + shape=[d1, d2 - len(ind)]) + return A diff --git a/caiman/utils/utils.py b/caiman/utils/utils.py index 4ea2f569b..9dca55cc3 100644 --- a/caiman/utils/utils.py +++ b/caiman/utils/utils.py @@ -343,11 +343,12 @@ def apply_magic_wand(A, gSig, dims, A_thr=None, coms=None, dview=None, return masks + def cell_magic_wand_wrapper(params): - a, com, min_radius, max_radius, roughness, zoom_factor, center_range = params - msk = cell_magic_wand(a, com, min_radius, max_radius, roughness, - zoom_factor, center_range) - return msk + a, com, min_radius, max_radius, roughness, zoom_factor, center_range = params + msk = cell_magic_wand(a, com, min_radius, max_radius, roughness, + zoom_factor, center_range) + return msk #%% From https://codereview.stackexchange.com/questions/120802/recursively-save-python-dictionaries-to-hdf5-files-using-h5py @@ -494,33 +495,3 @@ def recursively_load_dict_contents_from_group(h5file:h5py.File, path:str) -> Dic else: ans[key] = recursively_load_dict_contents_from_group(h5file, path + key + '/') return ans - - -def csc_column_remove(A: scipy.sparse.csc_matrix, ind: Iterable[int]) -> scipy.sparse.csc_matrix: - """ Removes specified columns for a scipy.sparse csc_matrix - Args: - A: scipy.sparse.csc_matrix - Input matrix - ind: iterable[int] - list or np.array with columns to be removed - """ - d1, d2 = A.shape - if 'csc_matrix' not in str(type(A)): - logging.warning("Original matrix not in csc_format. Converting it" + - " anyway.") - A = scipy.sparse.csc_matrix(A) - indptr = A.indptr - ind_diff = np.diff(A.indptr).tolist() - ind_sort = sorted(ind, reverse=True) - data_list = [A.data[indptr[i]:indptr[i+1]] for i in range(d2)] - indices_list = [A.indices[indptr[i]:indptr[i+1]] for i in range(d2)] - for i in ind_sort: - del data_list[i] - del indices_list[i] - del ind_diff[i] - indptr_final = np.cumsum([0] + ind_diff) - data_final = [item for sublist in data_list for item in sublist] - indices_final = [item for sublist in indices_list for item in sublist] - A = scipy.sparse.csc_matrix((data_final, indices_final, indptr_final), - shape=[d1, d2 - len(ind)]) - return A From 9df2433747332693ccd0f06df9e31464437d5a61 Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 24 Apr 2019 16:50:59 -0400 Subject: [PATCH 21/33] convert A to np.array prior to gt comparison --- caiman/source_extraction/cnmf/estimates.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index 949cfe152..5ef8bc65c 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -1066,7 +1066,7 @@ def remove_small_large_neurons(self, min_size_neuro, max_size_neuro): if self.A_thr is None: raise Exception('You need to compute thresolded components before calling remove_duplicates: use the threshold_components method') - A_gt_thr_bin = self.A_thr > 0 + A_gt_thr_bin = self.A_thr.toarray() > 0 size_neurons_gt = A_gt_thr_bin.sum(0) neurons_to_keep = np.where((size_neurons_gt > min_size_neuro) & (size_neurons_gt < max_size_neuro))[0] self.select_components(idx_components=neurons_to_keep) @@ -1088,7 +1088,7 @@ def remove_duplicates(self, predictions=None, r_values=None, dist_thr=0.1, min_d if self.A_thr is None: raise Exception('You need to compute thresolded components before calling remove_duplicates: use the threshold_components method') - A_gt_thr_bin = (self.A_thr > 0).reshape([self.dims[0], self.dims[1], -1], order='F').transpose([2, 0, 1]) * 1. + A_gt_thr_bin = (self.A_thr.toarray() > 0).reshape([self.dims[0], self.dims[1], -1], order='F').transpose([2, 0, 1]) * 1. duplicates_gt, indeces_keep_gt, indeces_remove_gt, D_gt, overlap_gt = detect_duplicates_and_subsets( A_gt_thr_bin,predictions=predictions, r_values=r_values,dist_thr=dist_thr, min_dist=min_dist, @@ -1138,8 +1138,8 @@ def compare_components(estimate_gt, estimate_cmp, Cn=None, thresh_cost=.8, min_ plt.figure(figsize=(20, 10)) dims = estimate_gt.dims - A_gt_thr_bin = (estimate_gt.A_thr>0).reshape([dims[0], dims[1], -1], order='F').transpose([2, 0, 1]) * 1. - A_thr_bin = (estimate_cmp.A_thr>0).reshape([dims[0], dims[1], -1], order='F').transpose([2, 0, 1]) * 1. + A_gt_thr_bin = (estimate_gt.A_thr.toarray()>0).reshape([dims[0], dims[1], -1], order='F').transpose([2, 0, 1]) * 1. + A_thr_bin = (estimate_cmp.A_thr.toarray()>0).reshape([dims[0], dims[1], -1], order='F').transpose([2, 0, 1]) * 1. tp_gt, tp_comp, fn_gt, fp_comp, performance_cons_off = nf_match_neurons_in_binary_masks( A_gt_thr_bin, A_thr_bin, thresh_cost=thresh_cost, min_dist=min_dist, print_assignment=print_assignment, From 682da148fda96108f5ddc7300ee46b4b0aa52ef3 Mon Sep 17 00:00:00 2001 From: barryza1 Date: Thu, 25 Apr 2019 11:45:21 -0400 Subject: [PATCH 22/33] add ignore preexisting cluster option --- caiman/cluster.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/caiman/cluster.py b/caiman/cluster.py index 6c6a22ff2..fb843e10c 100644 --- a/caiman/cluster.py +++ b/caiman/cluster.py @@ -74,7 +74,7 @@ def extract_patch_coordinates(dims, rf, stride, border_pix=0, indeces=[slice(Non sl_stop = [dim if sl.stop is None else sl.stop for (sl, dim) in zip(indeces, dims)] sl_step = [1 for sl in indeces] # not used dims_large = dims - dims = np.minimum(np.array(dims) - border_pix, sl_stop) - np.maximum(border_pix, sl_start) + dims = np.minimum(np.array(dims) - border_pix, sl_stop) - np.maximum(border_pix, sl_start) coords_flat = [] shapes = [] @@ -359,13 +359,16 @@ def stop_server(ipcluster='ipcluster', pdir=None, profile=None, dview=None): #%% -def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=False): +def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=False, ignore_preexisting=False): """Setup and/or restart a parallel cluster. Args: backend: str 'multiprocessing' [alias 'local'], 'ipyparallel', and 'SLURM' ipyparallel and SLURM backends try to restart if cluster running. backend='multiprocessing' raises an exception if a cluster is running. + ignore_preexisting: bool + If True, skips check for the existence of an already running multiprocessing + pool, which is usually indicative of a previously-started CaImAn cluster Returns: c: ipyparallel.Client object; only used for ipyparallel and SLURM backends, else None @@ -407,8 +410,13 @@ def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=Fal elif (backend == 'multiprocessing') or (backend == 'local'): if len(multiprocessing.active_children()) > 0: - raise Exception( - 'A cluster is already runnning. Terminate with dview.terminate() if you want to restart.') + if ignore_preexisting: + logger.warn('Found an existing multiprocessing pool. ' + 'This is often indicative of an already-running CaImAn cluster. ' + 'You have configured the cluster setup to not raise an exception.') + else: + raise Exception( + 'A cluster is already runnning. Terminate with dview.terminate() if you want to restart.') if (platform.system() == 'Darwin') and (sys.version_info > (3, 0)): try: if 'kernel' in get_ipython().trait_names(): # type: ignore @@ -423,7 +431,7 @@ def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=Fal except: # If we're not running under ipython, don't do anything. pass c = None - + dview = Pool(n_processes) else: raise Exception('Unknown Backend') From 96363c5add3a0d139d351d08e968a46f9e1f0fe7 Mon Sep 17 00:00:00 2001 From: barryza1 Date: Thu, 25 Apr 2019 11:46:45 -0400 Subject: [PATCH 23/33] update docstring --- caiman/cluster.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/caiman/cluster.py b/caiman/cluster.py index fb843e10c..f8e979db8 100644 --- a/caiman/cluster.py +++ b/caiman/cluster.py @@ -367,7 +367,7 @@ def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=Fal ipyparallel and SLURM backends try to restart if cluster running. backend='multiprocessing' raises an exception if a cluster is running. ignore_preexisting: bool - If True, skips check for the existence of an already running multiprocessing + If True, ignores the existence of an already running multiprocessing pool, which is usually indicative of a previously-started CaImAn cluster Returns: From 13939df70ccdfe537910fd5100541071793d6077 Mon Sep 17 00:00:00 2001 From: barryza1 Date: Sun, 28 Apr 2019 11:01:39 -0400 Subject: [PATCH 24/33] add repr to CNMFParams --- caiman/source_extraction/cnmf/params.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/caiman/source_extraction/cnmf/params.py b/caiman/source_extraction/cnmf/params.py index d7a50e58a..f065dbce6 100644 --- a/caiman/source_extraction/cnmf/params.py +++ b/caiman/source_extraction/cnmf/params.py @@ -8,6 +8,8 @@ from ...paths import caiman_datadir from .utilities import dict_compare, get_file_size +from pprint import pformat + class CNMFParams(object): def __init__(self, fnames=None, dims=None, dxy=(1, 1), @@ -700,7 +702,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'use_dense': use_dense, # flag for representation and storing of A and b 'use_peak_max': use_peak_max, # flag for finding candidate centroids } - + self.motion = { 'border_nan': 'copy', # flag for allowing NaN in the boundaries 'gSig_filt': None, # size of kernel for high pass spatial filtering in 1p data @@ -721,7 +723,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'upsample_factor_grid': 4, # motion field upsampling factor during FFT shifts 'use_cuda': False # flag for using a GPU } - + self.change_params(params_dict) if self.data['dims'] is None and self.data['fnames'] is not None: self.data['dims'] = get_file_size(self.data['fnames'], var_name_hdf5=self.data['var_name_hdf5'])[0] @@ -755,7 +757,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), "in group spatial automatically to False.") self.set('spatial', {'update_background_components': False}) if method_init=='corr_pnr' and ring_size_factor is not None: - logging.warning("using CNMF-E's ringmodel for background hence setting key " + + logging.warning("using CNMF-E's ringmodel for background hence setting key " + "normalize_init in group init automatically to False.") self.set('init', {'normalize_init': False}) @@ -842,6 +844,15 @@ def to_dict(self): 'merging': self.merging, 'motion': self.motion } + def __repr__(self): + + formatted_outputs = [ + '{}:\n\n{}'.format(group_name, pformat(group_dict)) + for group_name, group_dict in self.to_dict().items() + ] + + return 'CNMFParams:\n\n' + '\n\n'.join(formatted_outputs) + def change_params(self, params_dict, verbose=False): for gr in list(self.__dict__.keys()): self.set(gr, params_dict, verbose=verbose) From a72b3827a1100b6d9c5eae0e167fdbebab4d1378 Mon Sep 17 00:00:00 2001 From: epnev Date: Tue, 30 Apr 2019 07:26:59 -0400 Subject: [PATCH 25/33] calculating detrended traces for 1p data without normalizing --- caiman/source_extraction/cnmf/estimates.py | 11 ++++-- caiman/source_extraction/cnmf/utilities.py | 39 ++++++++++++++++++---- 2 files changed, 41 insertions(+), 9 deletions(-) diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index 949cfe152..88083cb53 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -579,7 +579,8 @@ def compute_residuals(self, Yr): return self def detrend_df_f(self, quantileMin=8, frames_window=500, - flag_auto=True, use_fast=False, use_residuals=True): + flag_auto=True, use_fast=False, use_residuals=True, + detrend_only=False): """Computes DF/F normalized fluorescence for the extracted traces. See caiman.source.extraction.utilities.detrend_df_f for details @@ -600,6 +601,11 @@ def detrend_df_f(self, quantileMin=8, frames_window=500, use_residuals: bool flag for using non-deconvolved traces in DF/F calculation + detrend_only: bool (False) + flag for only subtracting baseline and not normalizing by it. + Used in 1p data processing where baseline fluorescence cannot + be determined. + Returns: self: CNMF object self.F_dff contains the DF/F normalized traces @@ -623,7 +629,8 @@ def detrend_df_f(self, quantileMin=8, frames_window=500, self.F_dff = detrend_df_f(self.A, self.b, self.C, self.f, self.YrA, quantileMin=quantileMin, frames_window=frames_window, - flag_auto=flag_auto, use_fast=use_fast) + flag_auto=flag_auto, use_fast=use_fast, + detrend_only=detrend_only) return self def normalize_components(self): diff --git a/caiman/source_extraction/cnmf/utilities.py b/caiman/source_extraction/cnmf/utilities.py index 7895b1af7..d52681412 100644 --- a/caiman/source_extraction/cnmf/utilities.py +++ b/caiman/source_extraction/cnmf/utilities.py @@ -293,7 +293,7 @@ def extract_DF_F(Yr, A, C, bl, quantileMin=8, frames_window=200, block_size=400, return C_df def detrend_df_f(A, b, C, f, YrA=None, quantileMin=8, frames_window=500, - flag_auto=True, use_fast=False): + flag_auto=True, use_fast=False, detrend_only=False): """ Compute DF/F signal without using the original data. In general much faster than extract_DF_F @@ -323,7 +323,12 @@ def detrend_df_f(A, b, C, f, YrA=None, quantileMin=8, frames_window=500, flag for determining quantile automatically use_fast: bool - flag for using approximate fast percentile filtering + flag for u´sing approximate fast percentile filtering + + detrend_only: bool (False) + flag for only subtracting baseline and not normalizing by it. + Used in 1p data processing where baseline fluorescence cannot be + determined. Returns: F_df: @@ -333,7 +338,14 @@ def detrend_df_f(A, b, C, f, YrA=None, quantileMin=8, frames_window=500, if C is None: logging.warning("There are no components for DF/F extraction!") return None - + + if b is None or f is None: + b = np.zeros((A.shape[0], 1)) + f = np.zeros((1, C.shape[1])) + logging.warning("Background components not present. Results should" + + " not be interpreted as DF/F normalized but only" + + " as detrended.") + detrend_only = True if 'csc_matrix' not in str(type(A)): A = scipy.sparse.csc_matrix(A) if 'array' not in str(type(b)): @@ -362,7 +374,10 @@ def detrend_df_f(A, b, C, f, YrA=None, quantileMin=8, frames_window=500, zip(F, data_prct)]) Df = np.stack([np.percentile(f, prctileMin) for f, prctileMin in zip(B, data_prct)]) - F_df = (F - Fd[:, None]) / (Df[:, None] + Fd[:, None]) + if not detrend_only: + F_df = (F - Fd) / (Df[:, None] + Fd[:, None]) + else: + F_df = F - Fd else: if use_fast: Fd = np.stack([fast_prct_filt(f, level=prctileMin, @@ -378,18 +393,28 @@ def detrend_df_f(A, b, C, f, YrA=None, quantileMin=8, frames_window=500, Df = np.stack([scipy.ndimage.percentile_filter( f, prctileMin, (frames_window)) for f, prctileMin in zip(B, data_prct)]) - F_df = (F - Fd) / (Df + Fd) + if not detrend_only: + F_df = (F - Fd) / (Df + Fd) + else: + F_df = F - Fd else: if frames_window is None or frames_window > T: Fd = np.percentile(F, quantileMin, axis=1) Df = np.percentile(B, quantileMin, axis=1) - F_df = (F - Fd) / (Df[:, None] + Fd[:, None]) + if not detrend_only: + F_df = (F - Fd) / (Df[:, None] + Fd[:, None]) + else: + F_df = F - Fd else: Fd = scipy.ndimage.percentile_filter( F, quantileMin, (frames_window, 1)) Df = scipy.ndimage.percentile_filter( B, quantileMin, (frames_window, 1)) - F_df = (F - Fd) / (Df + Fd) + if not detrend_only: + F_df = (F - Fd) / (Df + Fd) + else: + F_df = F - Fd + return F_df def fast_prct_filt(input_data, level=8, frames_window=1000): From b862011f1ad95995e8419ba4fd6ce56a2cbc7afc Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Tue, 30 Apr 2019 12:22:51 -0400 Subject: [PATCH 26/33] Windows install docs: more flexible instructions on VS linkage --- INSTALL-windows.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/INSTALL-windows.md b/INSTALL-windows.md index 95745d9d3..e14533d2b 100755 --- a/INSTALL-windows.md +++ b/INSTALL-windows.md @@ -18,7 +18,7 @@ From that prompt. issue the following commands (if you wish to use the dev branc conda install -n caiman vs2017_win-64 ``` -At this point you will want to remove a startup script that visual studio made for your conda environment that can cause conda to crash while entering the caiman environment. Use the Windows find-file utility (under the Start Menu) to look for vs2015_compiler_vars.bat and/or vs2015_compiler_vars.bat under your home directory. At least one copy should show up. Delete the version that has conda\envs\caiman as part of its location. You may also want to do a search for keras_activate.bat under your home directory, find the one in conda\envs\caiman, and edit it so KERAS_BACKEND is set to tensorflow rather than theano. You may then continue the installation. +At this point you will want to remove a startup script that visual studio made for your conda environment that can cause conda to crash while entering the caiman environment. Use the Windows find-file utility (under the Start Menu) to look for vs2015_compiler_vars.bat and/or vs2015_compiler_vars.bat under your home directory. If the vs2015 prefix is not present, vs2018 or another year may. If you find one, delete the version that has conda\envs\caiman as part of its location. You may also want to do a search for keras_activate.bat under your home directory, find the one in conda\envs\caiman, and edit it so KERAS_BACKEND is set to tensorflow rather than theano. You may then continue the installation. ```bash activate caiman From a2fbb4e91495eb6a8de1f729fa39b8299e4469ac Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Tue, 30 Apr 2019 12:29:15 -0400 Subject: [PATCH 27/33] Windows install instructions: Clarify --- INSTALL-windows.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/INSTALL-windows.md b/INSTALL-windows.md index e14533d2b..7751a53c2 100755 --- a/INSTALL-windows.md +++ b/INSTALL-windows.md @@ -18,7 +18,7 @@ From that prompt. issue the following commands (if you wish to use the dev branc conda install -n caiman vs2017_win-64 ``` -At this point you will want to remove a startup script that visual studio made for your conda environment that can cause conda to crash while entering the caiman environment. Use the Windows find-file utility (under the Start Menu) to look for vs2015_compiler_vars.bat and/or vs2015_compiler_vars.bat under your home directory. If the vs2015 prefix is not present, vs2018 or another year may. If you find one, delete the version that has conda\envs\caiman as part of its location. You may also want to do a search for keras_activate.bat under your home directory, find the one in conda\envs\caiman, and edit it so KERAS_BACKEND is set to tensorflow rather than theano. You may then continue the installation. +At this point you will want to remove a startup script that visual studio made for your conda environment that can cause conda to crash while entering the caiman environment. Use the Windows find-file utility (under the Start Menu) to look for vs2015_compiler_vars.bat and/or vs2017_compiler_vars.bat under your home directory. If you find one, delete the version that has conda\envs\caiman as part of its location. Some systems versions will not have this file (this is fine and leaves you with nothing to do for this step). You may also want to do a search for keras_activate.bat under your home directory, find the one in conda\envs\caiman, and edit it so KERAS_BACKEND is set to tensorflow rather than theano. You may then continue the installation. ```bash activate caiman From 4d4d6df02521806dcf6d02deee706ac01ccbd7d1 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Tue, 30 Apr 2019 14:13:10 -0400 Subject: [PATCH 28/33] Adding unused.py as a "function purgatory" --- caiman/utils/unused.py | 70 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 caiman/utils/unused.py diff --git a/caiman/utils/unused.py b/caiman/utils/unused.py new file mode 100644 index 000000000..ae45773e2 --- /dev/null +++ b/caiman/utils/unused.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +####################### +# unused functions +# +# These are not presently used in the Caiman codebase, but are +# being kept for being potentially useful (and either not trivial to +# reimplement, or themselves a nontrivial idea). +# +# Code in this file should: +# 1) Not be considered for purposes of "don't break APIs in use". The +# person who would re-integrate the code has the responsibility for +# adapting it to current conventions +# 2) Never be called without moving it back out. This file should never be +# imported. +# 3) Be considered one step from deletion, if someone decides that the git +# history is good enough. +# 4) Note which file it came from +# +# Please prefer to delete code unless there's some reason to keep it + +# From caiman/cluster.py +def get_patches_from_image(img, shapes, overlaps): + # todo todocument + d1, d2 = np.shape(img) + rf = np.divide(shapes, 2) + _, coords_2d = extract_patch_coordinates(d1, d2, rf=rf, stride=overlaps) + imgs = np.empty(coords_2d.shape[:2], dtype=np.object) + + for idx_0, count_0 in enumerate(coords_2d): + for idx_1, count_1 in enumerate(count_0): + imgs[idx_0, idx_1] = img[count_1[0], count_1[1]] + + return imgs, coords_2d + +# From caiman/components_evaluation.py +def estimate_noise_mode(traces, robust_std=False, use_mode_fast=False, return_all=False): + """ estimate the noise in the traces under assumption that signals are sparse and only positive. The last dimension should be time. + + """ + # todo todocument + if use_mode_fast: + md = mode_robust_fast(traces, axis=1) + else: + md = mode_robust(traces, axis=1) + + ff1 = traces - md[:, None] + # only consider values under the mode to determine the noise standard deviation + ff1 = -ff1 * (ff1 < 0) + if robust_std: + # compute 25 percentile + ff1 = np.sort(ff1, axis=1) + ff1[ff1 == 0] = np.nan + Ns = np.round(np.sum(ff1 > 0, 1) * .5) + iqr_h = np.zeros(traces.shape[0]) + + for idx, _ in enumerate(ff1): + iqr_h[idx] = ff1[idx, -Ns[idx]] + + # approximate standard deviation as iqr/1.349 + sd_r = 2 * iqr_h / 1.349 + + else: + Ns = np.sum(ff1 > 0, 1) + sd_r = np.sqrt(old_div(np.sum(ff1**2, 1), Ns)) + + if return_all: + return md, sd_r + else: + return sd_r From 0b3110aee00c0402ba6326b8e4ad93a9b5e36361 Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Tue, 30 Apr 2019 14:24:55 -0400 Subject: [PATCH 29/33] Fix some "indeces" that slipped back in there --- caiman/source_extraction/cnmf/cnmf.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 1fbe9fcf1..b11b88af5 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -298,7 +298,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p self.estimates = Estimates(A=Ain, C=Cin, b=b_in, f=f_in, dims=self.params.data['dims']) - def fit_file(self, motion_correct=False, indeces=[slice(None)]*2): + def fit_file(self, motion_correct=False, indices=[slice(None)]*2): """ This method packages the analysis pipeline (motion correction, memory mapping, patch based CNMF processing) in a single method that can be @@ -311,7 +311,7 @@ def fit_file(self, motion_correct=False, indeces=[slice(None)]*2): Args: motion_correct (bool) flag for performing motion correction - indeces (list of slice objects) + indices (list of slice objects) perform analysis only on a part of the FOV Returns: cnmf object with the current estimates @@ -403,10 +403,10 @@ def fit(self, images, indices=[slice(None), slice(None)]): if len(indices) < len(images.shape): indices = indices + [slice(None)]*(len(images.shape) - len(indices)) dims_orig = images.shape[1:] - dims_sliced = images[tuple(indeces)].shape[1:] + dims_sliced = images[tuple(indices)].shape[1:] is_sliced = (dims_orig != dims_sliced) if self.params.get('patch', 'rf') is None and (is_sliced or 'ndarray' in str(type(images))): - images = images[tuple(indeces)] + images = images[tuple(indices)] self.dview = None logging.warning("Parallel processing in a single patch " "is not available for loaded in memory or sliced" + From 026909a57a592d143b150e30fce290798903dcff Mon Sep 17 00:00:00 2001 From: Pat Gunn Date: Tue, 30 Apr 2019 18:17:30 -0400 Subject: [PATCH 30/33] Fix comparison tests. Add docs on managing tensorflow/opencv versions on dated distros. --- README-Distros.md | 27 ++++++++++++++++++++++++ README.md | 1 + caiman/tests/comparison_humans.py | 2 +- caiman/tests/comparison_humans_online.py | 2 +- 4 files changed, 30 insertions(+), 2 deletions(-) create mode 100644 README-Distros.md diff --git a/README-Distros.md b/README-Distros.md new file mode 100644 index 000000000..eb7972e0e --- /dev/null +++ b/README-Distros.md @@ -0,0 +1,27 @@ +CaImAn and Older Linux Distros +===================== +On some older Linux distros, you may have trouble building a working environment because some Conda packages are introducing glibc version requirements into their binary builds. This has been seen on RHEL7/CentOS7 and may show up in other distributions as well. The errors turn up when tensorflow is imported from python code, with the loader failing because it cannot find symbols indicating a modern libc. + +Fixing this is not trivial, but it is doable. + +Fixing it +========= +To resolve this, you will first want to install a different build of Tensorflow, not from conda-forge. I used +``` +conda search -c conda-forge tensorflow +``` + +to find a version from the main distro built against mkl (you will see mkl in the build string). Having done this, I installed that particular build, in my case by doing the following (your steps may vary): +``` +conda install tensorflow=1.13.1=mkl_py37h54b294f_0 +``` + +This works, but in doing so it changed your version of opencv and a number of other things to non-conda-forge versions. The non-conda-forge builds of opencv are built without graphical bindings, making them not useful for some things in CaImAn. We can switch them back by looking for available builds of opencv using a conda search command like the above, and then selecting a (similar version to what we got) suitable conda-forge build of opencv. In my case I finished fixing it with: + +``` +conda install -c conda-forge opencv=3.4.4 +``` + +Support +======= +This is not a pretty procedure, but we will try to support it if it is necessary in your environment. diff --git a/README.md b/README.md index abc5ecbaa..5d2bd3eb9 100755 --- a/README.md +++ b/README.md @@ -117,6 +117,7 @@ If you prefer to manage this information somewhere else, the `CAIMAN_DATA` envir Alternative environments: * [Using GPU](/README-GPU.md) + * [Older Linux Distros](/README-Distros.md) ### Known Issues diff --git a/caiman/tests/comparison_humans.py b/caiman/tests/comparison_humans.py index 229909a5d..1eae69bf6 100644 --- a/caiman/tests/comparison_humans.py +++ b/caiman/tests/comparison_humans.py @@ -400,7 +400,7 @@ cnm2.estimates.plot_contours(img=Cn) # %% prepare ground truth masks gt_file = os.path.join(os.path.split(fname_new)[0], os.path.split(fname_new)[1][:-4] + 'match_masks.npz') - with np.load(gt_file, encoding='latin1') as ld: + with np.load(gt_file, encoding='latin1', allow_pickle=True) as ld: print(ld.keys()) Cn_orig = ld['Cn'] diff --git a/caiman/tests/comparison_humans_online.py b/caiman/tests/comparison_humans_online.py index a26e76180..16aeb7714 100644 --- a/caiman/tests/comparison_humans_online.py +++ b/caiman/tests/comparison_humans_online.py @@ -244,7 +244,7 @@ # %% load consensus annotations and filter for size gt_file = glob.glob(os.path.join(base_folder, params_movie[ind_dataset]['folder_name'], '*masks.npz'))[0] - with np.load(gt_file, encoding='latin1') as ld: + with np.load(gt_file, encoding='latin1', allow_pickle=True) as ld: d1_or = int(ld['d1']) d2_or = int(ld['d2']) dims_or = (d1_or, d2_or) From 0343a5176c2bd44f283201ed6d897a0411b7f1d8 Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 1 May 2019 11:21:15 -0400 Subject: [PATCH 31/33] add allow_pickle flag in a bunch of places --- use_cases/eLife_scripts/Figure_4-figure_supplement1.py | 2 +- use_cases/eLife_scripts/figure_4_a_c_d_e.py | 4 ++-- use_cases/eLife_scripts/figure_4_b.py | 4 ++-- use_cases/eLife_scripts/figure_5.py | 6 +++--- .../preprocessing_files/Preprocess_CaImAn_batch.py | 6 +++--- .../preprocessing_files/Preprocess_CaImAn_online.py | 2 +- 6 files changed, 12 insertions(+), 12 deletions(-) diff --git a/use_cases/eLife_scripts/Figure_4-figure_supplement1.py b/use_cases/eLife_scripts/Figure_4-figure_supplement1.py index 6e06b4ea3..e294577de 100644 --- a/use_cases/eLife_scripts/Figure_4-figure_supplement1.py +++ b/use_cases/eLife_scripts/Figure_4-figure_supplement1.py @@ -20,7 +20,7 @@ base_folder = '/mnt/ceph/neuro/DataForPublications/DATA_PAPER_ELIFE/WEBSITE' -with np.load(os.path.join(base_folder, 'all_records_grid_online.npz')) as ld: +with np.load(os.path.join(base_folder, 'all_records_grid_online.npz'), allow_pickle=True) as ld: records_online = ld['records'] records_online = [list(rec) for rec in records_online] diff --git a/use_cases/eLife_scripts/figure_4_a_c_d_e.py b/use_cases/eLife_scripts/figure_4_a_c_d_e.py index 90e6c30e4..ef184c980 100644 --- a/use_cases/eLife_scripts/figure_4_a_c_d_e.py +++ b/use_cases/eLife_scripts/figure_4_a_c_d_e.py @@ -96,9 +96,9 @@ def precision_snr(snr_gt, snr_gt_fn, snr_cnmf, snr_cnmf_fp, snr_thrs): # %% RELOAD ALL THE RESULTS INSTEAD OF REGENERATING THEM -with np.load(os.path.join(base_folder, 'all_res_web.npz')) as ld: +with np.load(os.path.join(base_folder, 'all_res_web.npz'), allow_pickle=True) as ld: all_results = ld['all_results'][()] -with np.load(os.path.join(base_folder, 'all_res_online_web_bk.npz')) as ld: +with np.load(os.path.join(base_folder, 'all_res_online_web_bk.npz'), allow_pickle=True) as ld: all_results_online = ld['all_results'][()] pl.rcParams['pdf.fonttype'] = 42 diff --git a/use_cases/eLife_scripts/figure_4_b.py b/use_cases/eLife_scripts/figure_4_b.py index e27900b5e..c0a58645b 100644 --- a/use_cases/eLife_scripts/figure_4_b.py +++ b/use_cases/eLife_scripts/figure_4_b.py @@ -34,7 +34,7 @@ base_folder = '/mnt/ceph/neuro/DataForPublications/DATA_PAPER_ELIFE/WEBSITE/' # %% Figure 4b and GRID statistics -with np.load(os.path.join(base_folder,'ALL_RECORDS_GRID_FINAL.npz')) as ld: +with np.load(os.path.join(base_folder,'ALL_RECORDS_GRID_FINAL.npz'), allow_pickle=True) as ld: records = ld['records'][()] records = [list(rec) for rec in records] records = [rec[:5]+[float(rr) for rr in rec[5:]] for rec in records] @@ -42,7 +42,7 @@ columns = ['name', 'gr_snr', 'grid_rval', 'grid_max_prob_rej', 'grid_thresh_CNN', 'recall', 'precision', 'f1_score'] -with np.load(os.path.join(base_folder, 'all_records_grid_online.npz')) as ld: +with np.load(os.path.join(base_folder, 'all_records_grid_online.npz'), allow_pickle=True) as ld: records_online = ld['records'] records_online = [list(rec) for rec in records_online] records_online = [rec[:4] + [float(rr) for rr in rec[4:]] for rec in records_online] diff --git a/use_cases/eLife_scripts/figure_5.py b/use_cases/eLife_scripts/figure_5.py index 9ae033360..9adb99c66 100644 --- a/use_cases/eLife_scripts/figure_5.py +++ b/use_cases/eLife_scripts/figure_5.py @@ -50,15 +50,15 @@ fls = glob.glob('.') # here should go the path to the generated files all_results_online = dict() for fl in fls: - with np.load(fl) as ld: + with np.load(fl, allow_pickle=True) as ld: all_results_online.update((ld['all_results'][()])) np.savez(os.path.join(base_folder, 'all_res_online_web.npz'), all_results=all_results_online) # %% RELOAD ALL THE RESULTS INSTEAD OF REGENERATING THEM -with np.load(os.path.join(base_folder, 'all_res_web.npz')) as ld: +with np.load(os.path.join(base_folder, 'all_res_web.npz'), allow_pickle=True) as ld: all_results = ld['all_results'][()] -with np.load(os.path.join(base_folder, 'all_res_online_web_bk.npz')) as ld: +with np.load(os.path.join(base_folder, 'all_res_online_web_bk.npz'), allow_pickle=True) as ld: all_results_online = ld['all_results'][()] pl.rcParams['pdf.fonttype'] = 42 diff --git a/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_batch.py b/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_batch.py index baf6b7bac..e6b76aa84 100644 --- a/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_batch.py +++ b/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_batch.py @@ -407,7 +407,7 @@ cnm2.estimates.plot_contours(img=Cn) # %% prepare ground truth masks gt_file = os.path.join(os.path.split(fname_new)[0], os.path.split(fname_new)[1][:-4] + 'match_masks.npz') - with np.load(gt_file, encoding='latin1') as ld: + with np.load(gt_file, encoding='latin1', allow_pickle=True) as ld: print(ld.keys()) Cn_orig = ld['Cn'] @@ -476,7 +476,7 @@ performance_tmp['t_refit'] = t_refit performance_tmp['t_eva'] = t_eva_comps - with np.load(os.path.join(base_folder,fname_new.split('/')[-2],'gt_eval.npz')) as ld: + with np.load(os.path.join(base_folder,fname_new.split('/')[-2],'gt_eval.npz'), allow_pickle=True) as ld: print(ld.keys()) performance_tmp.update(ld) @@ -508,7 +508,7 @@ npzfile = os.path.join(base_folder,params_movie['fname'][:-5] + '_perf_web_gsig.npz') print(npzfile) if os.path.exists(npzfile): - with np.load(npzfile) as ld: + with np.load(npzfile, allow_pickle=True) as ld: all_results[params_movie['fname'].split('/')[-2]] = ld['all_results'][()] else: print("*** NOT EXIST ***" + npzfile) diff --git a/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_online.py b/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_online.py index e570b52bd..4883c7a88 100644 --- a/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_online.py +++ b/use_cases/eLife_scripts/preprocessing_files/Preprocess_CaImAn_online.py @@ -286,7 +286,7 @@ # %% load consensus annotations and filter for size gt_file = glob.glob(os.path.join(base_folder, params_movie[ind_dataset]['folder_name'], '*masks.npz'))[0] - with np.load(gt_file, encoding='latin1') as ld: + with np.load(gt_file, encoding='latin1', allow_pickle=True) as ld: d1_or = int(ld['d1']) d2_or = int(ld['d2']) dims_or = (d1_or, d2_or) From 86335f90f65801916d08c15539c0dfa99e0b15f0 Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 1 May 2019 11:48:15 -0400 Subject: [PATCH 32/33] roll back to 3.6 --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 013a1eb0e..338e37249 100644 --- a/environment.yml +++ b/environment.yml @@ -1,7 +1,7 @@ channels: - conda-forge dependencies: -- python=3.7 +- python=3.6 - tensorflow # Must list before keras to preempt default keras backend - bokeh - cython From 918f82cbb9ac68c285de74a38722167d43483e0a Mon Sep 17 00:00:00 2001 From: epnev Date: Wed, 1 May 2019 13:39:56 -0400 Subject: [PATCH 33/33] add message when removing empty components --- caiman/source_extraction/cnmf/temporal.py | 1 + 1 file changed, 1 insertion(+) diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 28e20479e..29d0c4f4a 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -216,6 +216,7 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N ff = np.where(np.sum(C, axis=1) == 0) # remove empty components if np.size(ff) > 0: # Eliminating empty temporal components ff = ff[0] + logging.info('removing {0} empty spatial component(s)'.format(len(ff))) keep = list(range(A.shape[1])) for i in ff: keep.remove(i)