diff --git a/INSTALL-windows.md b/INSTALL-windows.md index 95745d9d3..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. 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 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 diff --git a/Jenkinsfile b/Jenkinsfile index 8524e011b..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\\vs2015_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' } } 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 c25ce9959..babc44ff4 100755 --- a/README.md +++ b/README.md @@ -123,6 +123,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/base/movies.py b/caiman/base/movies.py index a03b2d5cb..bb7a24831 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, Any]: """ 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) -> Any: """ - 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: @@ -1374,10 +1397,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) -> Any: """ load movies from list of file names Args: @@ -1393,7 +1416,7 @@ def load_movie_chain(file_list, fr=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 """ @@ -1422,47 +1445,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 @@ -1510,11 +1533,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: @@ -1564,11 +1589,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: @@ -1596,20 +1622,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) -> Any: ''' + 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: + movie ''' mov:List = [] print('unzipping file into movie object') @@ -1619,7 +1646,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)) @@ -1636,13 +1663,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/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 6c6a22ff2..ba83c5fa6 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, indices=[slice(None)]*2) -> Tuple[List, List]: """ Partition the FOV in patches and return the indexed in 2D and 1D (flatten, order='F') formats @@ -70,11 +56,11 @@ def extract_patch_coordinates(dims, rf, stride, border_pix=0, indeces=[slice(Non 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) + dims = np.minimum(np.array(dims) - border_pix, sl_stop) - np.maximum(border_pix, sl_start) coords_flat = [] shapes = [] @@ -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,13 +351,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:str='multiprocessing', n_processes:int=None, single_thread:bool=False, ignore_preexisting:bool=False) -> Tuple[Any, Any, Optional[int]]: """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, ignores 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 +402,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 +423,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') diff --git a/caiman/components_evaluation.py b/caiman/components_evaluation.py index 73bd8fc23..b6a3cbab7 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[np.array, np.array, 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) @@ -529,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): + 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 @@ -582,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=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) -> 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. @@ -643,6 +591,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 5de9ac792..a4dade990 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)') @@ -375,8 +378,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 +390,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) diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 9b1ac9e20..b11b88af5 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 @@ -297,7 +298,24 @@ 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 + 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 + indices (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] @@ -333,7 +351,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): """ @@ -355,17 +373,16 @@ 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. - - 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: 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 @@ -380,16 +397,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[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[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" + @@ -509,7 +526,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() @@ -540,7 +557,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") @@ -569,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, @@ -581,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) @@ -606,11 +622,12 @@ 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 - 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,\ @@ -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 @@ -872,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, @@ -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/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index 949cfe152..8d4e85e93 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 @@ -80,10 +82,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 @@ -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 @@ -579,7 +581,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 +603,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 +631,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): @@ -669,10 +678,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 +788,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 +815,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 +892,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 +927,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 +1009,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 +1019,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], @@ -1066,7 +1075,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,9 +1097,9 @@ 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( + 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 +1110,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]) @@ -1138,8 +1147,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, diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index fe9b921c7..cc25cc72b 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -34,6 +34,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 .deconvolution import constrained_foopsi @@ -634,7 +635,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) @@ -645,29 +646,29 @@ 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() + 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] diff --git a/caiman/source_extraction/cnmf/map_reduce.py b/caiman/source_extraction/cnmf/map_reduce.py index 719e24500..e034b5bd8 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 @@ -144,7 +145,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 +227,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): @@ -427,7 +428,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): @@ -461,9 +466,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 5b162cd67..852c274de 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -132,14 +132,14 @@ 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: + 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) @@ -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)) @@ -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] @@ -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): diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index a043a3785..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 @@ -1466,34 +1467,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) @@ -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: @@ -1589,7 +1590,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 +1600,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 +1621,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 +1654,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 +1668,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 +1685,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 +1697,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 +1766,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 @@ -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') 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) diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index b3473eb74..9547c6cb4 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.stats import csc_column_remove def basis_denoising(y, c, boh, sn, id2_, px): @@ -208,10 +209,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]) - A_ = np.zeros((d, nr + np.size(f, 0))) # init A_ + 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))) if dview is not None: if 'multiprocessing' in str(type(dview)): parallel_result = dview.map_async( @@ -222,22 +224,31 @@ 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, 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) + 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: if low_rank_background: background_ff = list(filter(lambda i: i >= nb, ff - nr)) @@ -246,16 +257,17 @@ 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") 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: @@ -487,8 +499,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,12 +516,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 @@ -586,7 +602,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: @@ -604,7 +620,7 @@ def threshold_components_parallel(pars): BW = BW.flatten() Ath2[BW] = Ath[BW] - return Ath2, i + return csr_matrix(Ath2), i # %% @@ -843,7 +859,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.csc_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() @@ -889,6 +906,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]) @@ -900,7 +920,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...') @@ -918,8 +943,12 @@ 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 + 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') @@ -927,7 +956,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): """ diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 86a7610e2..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) @@ -389,11 +390,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 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): diff --git a/caiman/tests/comparison/comparison.py b/caiman/tests/comparison/comparison.py index c8382df23..5552aab14 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'][()] @@ -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 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) 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/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 diff --git a/caiman/utils/utils.py b/caiman/utils/utils.py index b7eb41ee9..9dca55cc3 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) @@ -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 @@ -493,4 +494,4 @@ 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 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 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)