diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml index 7a10415ad..2904dacc1 100644 --- a/.github/workflows/run_tests.yml +++ b/.github/workflows/run_tests.yml @@ -26,7 +26,7 @@ jobs: with: auto-update-conda: true python-version: ${{ matrix.python-version }} - environment-file: environment.yml + environment-file: environment-minimal.yml activate-environment: caiman - name: Install OS Dependencies diff --git a/README.md b/README.md index 7c850c993..36d57acc6 100644 --- a/README.md +++ b/README.md @@ -62,6 +62,8 @@ The main use cases and notebooks are listed in the following table: A comprehensive list of references, where you can find detailed discussion of the methods and their development, can be found [here](https://caiman.readthedocs.io/en/master/CaImAn_features_and_references.html#references). +# CLI demos +Caiman also provides commandline demos, similar to the notebooks, demonstrating how to work with the codebase outside of Jupyter. They take their configuration primarily from json files (which you will want to modify to work with your data and its specifics) and should be reasonably easy to modify if they don't already do what you want them to do (in particular, saving things; a standard output format for Caiman is something intended for future releases). To run them, activate your environment, and find the demos in demos/general under your caiman data directory; you can run them like you would any other python application, or edit them with your code editor. Each demo comes with a json configuration file that you can customise. There is a README in the demos directory that covers some of this. # How to get help - [Online documentation](https://caiman.readthedocs.io/en/latest/) contains a lot of general information about Caiman, the parameters, how to interpret its outputs, and more. diff --git a/caiman/base/movies.py b/caiman/base/movies.py index 60162a52e..9f44dac43 100644 --- a/caiman/base/movies.py +++ b/caiman/base/movies.py @@ -673,69 +673,6 @@ def NonnegativeMatrixFactorization(self, return space_components, time_components - 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. - Implements bith the nmf and nnsc methods - - Args: - n_components: int - - method: 'nnsc' or 'nmf' (see http://spams-devel.gforge.inria.fr/doc-python/html/index.html) - - lambda1: see http://spams-devel.gforge.inria.fr/doc-python/html/index.html - - iterations: see http://spams-devel.gforge.inria.fr/doc-python/html/index.html - - batchsize: see http://spams-devel.gforge.inria.fr/doc-python/html/index.html - - model: see http://spams-devel.gforge.inria.fr/doc-python/html/index.html - - **kwargs: more arguments to be passed to nmf or nnsc - - Returns: - time_comps - - space_comps - """ - try: - import spams # XXX consider moving this to the head of the file - except: - logging.error("You need to install the SPAMS package") - raise - - T, d1, d2 = np.shape(self) - d = d1 * d2 - X = np.asfortranarray(np.reshape(self, [T, d], order='F')) - - if method == 'nmf': - (time_comps, V) = spams.nmf(X, return_lasso=True, K=n_components, numThreads=4, iter=iterations, **kwargs) - - elif method == 'nnsc': - (time_comps, V) = spams.nnsc(X, - return_lasso=True, - K=n_components, - lambda1=lambda1, - iter=iterations, - model=model, - **kwargs) - else: - raise Exception('Method unknown') - - space_comps = [] - - for _, mm in enumerate(V): - space_comps.append(np.reshape(mm.todense(), (d1, d2), order='F')) - - return time_comps, np.array(space_comps) - 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 @@ -1244,7 +1181,8 @@ def load(file_name: Union[str, list[str]], dimension of the movie along x and y if loading from a two dimensional numpy array var_name_hdf5: str - if loading from hdf5/n5 name of the dataset inside the file to load (ignored if the file only has one dataset) + if loading from hdf5/n5 name of the dataset inside the file to load (ignored if the file only has one dataset). + This is also used for (new-style) mat files in_memory: bool=False This changes the behaviour of the function for npy files to be a readwrite rather than readonly memmap, @@ -1314,17 +1252,6 @@ def load(file_name: Union[str, list[str]], basename, extension = os.path.splitext(file_name) extension = extension.lower() - if extension == '.mat': - logging.warning('Loading a *.mat file. x- and y- dimensions ' + - 'might have been swapped.') - try: # scipy >= 1.8 - byte_stream, file_opened = scipy.io.matlab._mio._open_file(file_name, appendmat=False) - mjv, mnv = scipy.io.matlab.miobase.get_matfile_version(byte_stream) - except: # scipy <= 1.7 - byte_stream, file_opened = scipy.io.matlab.mio._open_file(file_name, appendmat=False) - mjv, mnv = scipy.io.matlab.mio.get_matfile_version(byte_stream) - if mjv == 2: - extension = '.h5' if extension in ['.tif', '.tiff', '.btf']: # load tif file with tifffile.TiffFile(file_name) as tffl: @@ -1512,23 +1439,23 @@ def load(file_name: Union[str, list[str]], else: input_arr = input_arr[np.newaxis, :, :] - elif extension == '.mat': # load npy file - input_arr = scipy.io.loadmat(file_name)['data'] - input_arr = np.rollaxis(input_arr, 2, -3) - if subindices is not None: - input_arr = input_arr[subindices] - elif extension == '.npz': # load movie from saved file if subindices is not None: raise Exception('Subindices not implemented') with np.load(file_name) as f: return movie(**f).astype(outtype) - elif extension in ('.hdf5', '.h5', '.nwb', 'n5', 'zarr'): + elif extension in ('.hdf5', '.h5', '.mat', '.nwb', 'n5', 'zarr'): if extension in ('n5', 'zarr'): # Thankfully, the zarr library lines up closely with h5py past the initial open f = zarr.open(file_name, "r") else: - f = h5py.File(file_name, "r") + try: + f = h5py.File(file_name, "r") + except: + if extension == '.mat': + raise Exception(f"Problem loading {file_name}: Unknown format. This may be in the original version 1 (non-hdf5) mat format; please convert it first") + else: + raise Exception(f"Problem loading {file_name}: Unknown format.") ignore_keys = ['__DATA_TYPES__'] # Known metadata that tools provide, add to this as needed. Sync with get_file_size() !! fkeys = list(filter(lambda x: x not in ignore_keys, f.keys())) if len(fkeys) == 1: # If the file we're parsing has only one dataset inside it, @@ -1951,11 +1878,17 @@ def load_iter(file_name: Union[str, list[str]], subindices=None, var_name_hdf5: yield frame # was frame[..., 0].astype(outtype) return - elif extension in ('.hdf5', '.h5', '.nwb', '.mat', 'n5', 'zarr'): - if extension in ('n5', 'zarr'): # Thankfully, the zarr library lines up closely with h5py past the initial open + elif extension in ('.hdf5', '.h5', '.nwb', '.mat', '.n5', '.zarr'): + if extension in ('.n5', '.zarr'): # Thankfully, the zarr library lines up closely with h5py past the initial open f = zarr.open(file_name, "r") else: - f = h5py.File(file_name, "r") + try: + f = h5py.File(file_name, "r") + except: + if extension == '.mat': + raise Exception(f"Problem loading {file_name}: Unknown format. This may be in the original version 1 (non-hdf5) mat format; please convert it first") + else: + raise Exception(f"Problem loading {file_name}: Unknown format.") ignore_keys = ['__DATA_TYPES__'] # Known metadata that tools provide, add to this as needed. fkeys = list(filter(lambda x: x not in ignore_keys, f.keys())) if len(fkeys) == 1: # If the hdf5 file we're parsing has only one dataset inside it, @@ -2010,11 +1943,7 @@ def get_file_size(file_name, var_name_hdf5:str='mov') -> tuple[tuple, Union[int, if os.path.exists(file_name): _, extension = os.path.splitext(file_name)[:2] extension = extension.lower() - if extension == '.mat': - byte_stream, file_opened = scipy.io.matlab.mio._open_file(file_name, appendmat=False) - mjv, mnv = scipy.io.matlab.mio.get_matfile_version(byte_stream) - if mjv == 2: - extension = '.h5' + if extension in ['.tif', '.tiff', '.btf']: tffl = tifffile.TiffFile(file_name) siz = tffl.series[0].shape @@ -2042,12 +1971,18 @@ def get_file_size(file_name, var_name_hdf5:str='mov') -> tuple[tuple, Union[int, filename = os.path.split(file_name)[-1] Yr, dims, T = caiman.mmapping.load_memmap(os.path.join( os.path.split(file_name)[0], filename)) - elif extension in ('.h5', '.hdf5', '.nwb', 'n5', 'zarr'): + elif extension in ('.h5', '.hdf5', '.mat', '.nwb', 'n5', 'zarr'): # FIXME this doesn't match the logic in load() if extension in ('n5', 'zarr'): # Thankfully, the zarr library lines up closely with h5py past the initial open f = zarr.open(file_name, "r") else: - f = h5py.File(file_name, "r") + try: + f = h5py.File(file_name, "r") + except: + if extension == '.mat': + raise Exception(f"Problem loading {file_name}: Unknown format. This may be in the original version 1 (non-hdf5) mat format; please convert it first") + else: + raise Exception(f"Problem loading {file_name}: Unknown format.") ignore_keys = ['__DATA_TYPES__'] # Known metadata that tools provide, add to this as needed. Sync with movies.my:load() !! kk = list(filter(lambda x: x not in ignore_keys, f.keys())) if len(kk) == 1: # TODO: Consider recursing into a group to find a dataset diff --git a/caiman/base/timeseries.py b/caiman/base/timeseries.py index f8868cb1c..55aa127d5 100644 --- a/caiman/base/timeseries.py +++ b/caiman/base/timeseries.py @@ -147,6 +147,7 @@ def save(self, Args: file_name: str name of file. Possible formats are tif, avi, npz, mmap and hdf5 + If a path is not part of the filename, it will be saved into a temporary directory under caiman_data to32: Bool whether to transform to 32 bits @@ -165,6 +166,9 @@ def save(self, if saving as .tif, specifies the compression level if saving as .avi or .mkv, compress=0 uses the IYUV codec, otherwise the FFV1 codec is used + Returns: + generated_filename: The full filename, path included, where the data was saved + Raises: Exception 'Extension Unknown' @@ -197,6 +201,8 @@ def foo(i): if to32 and not ('float32' in str(self.dtype)): curfr = curfr.astype(np.float32) tif.save(curfr, compress=compress) + return file_name + elif extension == '.npz': if to32 and not ('float32' in str(self.dtype)): input_arr = self.astype(np.float32) @@ -209,6 +215,8 @@ def foo(i): fr=self.fr, meta_data=self.meta_data, file_name=self.file_name) + return file_name + elif extension in ('.avi', '.mkv'): codec = None if compress == 0: @@ -241,6 +249,7 @@ def foo(i): for d in data: vw.write(cv2.cvtColor(d, cv2.COLOR_GRAY2BGR)) vw.release() + return file_name elif extension == '.mat': if self.file_name[0] is not None: @@ -271,6 +280,7 @@ def foo(i): 'meta_data': self.meta_data, 'file_name': f_name }) + return file_name elif extension in ('.hdf5', '.h5'): with h5py.File(file_name, "w") as f: @@ -289,6 +299,7 @@ def foo(i): if self.meta_data[0] is not None: logging.debug("Metadata for saved file: " + str(self.meta_data)) dset.attrs["meta_data"] = cpk.dumps(self.meta_data) + return file_name elif extension == '.mmap': base_name = name diff --git a/caiman/behavior/behavior.py b/caiman/behavior/behavior.py index ce9fc87be..9863f5603 100644 --- a/caiman/behavior/behavior.py +++ b/caiman/behavior/behavior.py @@ -37,6 +37,10 @@ def select_roi(img: np.ndarray, n_rois: int = 1) -> list: each element is an the mask considered a ROIs """ + # FIXME This function depends on particular builds of OpenCV + # and may be difficult to support moving forward; would be good to + # move this kind of code out of the core and find more portable ways + # to do it masks = [] for _ in range(n_rois): fig = plt.figure() @@ -130,8 +134,8 @@ def extract_magnitude_and_angle_from_OF(spatial_filter_, x, y = scipy.signal.medfilt(time_trace, kernel_size=[1, 1]).T x = scipy.signal.savgol_filter(x.squeeze(), sav_filter_size, 1) y = scipy.signal.savgol_filter(y.squeeze(), sav_filter_size, 1) - mag, dirct = to_polar(x - caiman.components_evaluation.mode_robust(x), - y - caiman.components_evaluation.mode_robust(y)) + mag, dirct = to_polar(x - caiman.utils.stats.mode_robust(x), + y - caiman.utils.stats.mode_robust(y)) dirct = scipy.signal.medfilt(dirct.squeeze(), kernel_size=1).T # normalize to pixel units @@ -325,25 +329,12 @@ def extract_components(mov_tot, if method_factorization == 'nmf': nmf = NMF(n_components=n_components, **kwargs) - time_trace = nmf.fit_transform(newm) spatial_filter = nmf.components_ spatial_filter = np.concatenate([np.reshape(sp, (d1, d2))[np.newaxis, :, :] for sp in spatial_filter], axis=0) - - elif method_factorization == 'dict_learn': - import spams - newm = np.asfortranarray(newm, dtype=np.float32) - time_trace = spams.trainDL(newm, K=n_components, mode=0, lambda1=1, posAlpha=True, iter=max_iter_DL) - - spatial_filter = spams.lasso(newm, - D=time_trace, - return_reg_path=False, - lambda1=0.01, - mode=spams.spams_wrap.PENALTY, - pos=True) - - spatial_filter = np.concatenate([np.reshape(sp, (d1, d2))[np.newaxis, :, :] for sp in spatial_filter.toarray()], - axis=0) + else: + # Caiman used to support a method_factorization called dict_learn, implemented using spams.lasso + raise Exception(f"Unknown or unsupported method_factorization: {method_factorization}") time_trace = [np.reshape(ttr, (c, T)).T for ttr in time_trace.T] diff --git a/caiman/caimanmanager.py b/caiman/caimanmanager.py index 457cb5082..c2d57b8c0 100755 --- a/caiman/caimanmanager.py +++ b/caiman/caimanmanager.py @@ -1,7 +1,6 @@ #!/usr/bin/env python import argparse -import distutils.dir_util import filecmp import glob import os @@ -53,20 +52,21 @@ def do_install_to(targdir: str, inplace: bool = False, force: bool = False) -> None: global sourcedir_base + ignore_pycache=shutil.ignore_patterns('__pycache__') if os.path.isdir(targdir) and not force: raise Exception(targdir + " already exists. You may move it out of the way, remove it, or use --force") if not inplace: # In this case we rely on what setup.py put in the share directory for the module if not force: - shutil.copytree(sourcedir_base, targdir) + shutil.copytree(sourcedir_base, targdir, ignore=ignore_pycache) else: - distutils.dir_util.copy_tree(sourcedir_base, targdir) + shutil.copytree(sourcedir_base, targdir, ignore=ignore_pycache, dirs_exist_ok=True) os.makedirs(os.path.join(targdir, 'temp' ), exist_ok=True) else: # here we recreate the other logical path here. Maintenance concern: Keep these reasonably in sync with what's in setup.py for copydir in extra_dirs: if not force: - shutil.copytree(copydir, os.path.join(targdir, copydir)) + shutil.copytree(copydir, os.path.join(targdir, copydir), ignore=ignore_pycache) else: - distutils.dir_util.copy_tree(copydir, os.path.join(targdir, copydir)) + shutil.copytree(copydir, os.path.join(targdir, copydir), ignore=ignore_pycache, dirs_exist_ok=True) os.makedirs(os.path.join(targdir, 'example_movies'), exist_ok=True) os.makedirs(os.path.join(targdir, 'temp' ), exist_ok=True) for stdmovie in standard_movies: diff --git a/caiman/mmapping.py b/caiman/mmapping.py index 65dd10cf4..ce4b040e8 100644 --- a/caiman/mmapping.py +++ b/caiman/mmapping.py @@ -405,7 +405,7 @@ def save_memmap(filenames:list[str], recompute_each_memmap = True - if recompute_each_memmap or (remove_init>0) or (idx_xy is not None)\ + if recompute_each_memmap or (remove_init > 0) or (idx_xy is not None)\ or (xy_shifts is not None) or (add_to_movie != 0) or (border_to_0>0)\ or slices is not None: @@ -527,7 +527,7 @@ def save_memmap(filenames:list[str], sys.stdout.flush() Ttot = Ttot + T - fname_new = caiman.paths.fn_relocated(fname_tot + f'_frames_{Ttot}.mmap') + fname_new = os.path.join(caiman.paths.get_tempdir(), caiman.paths.fn_relocated(f'{fname_tot}_frames_{Ttot}.mmap')) try: # need to explicitly remove destination on windows os.unlink(fname_new) diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py index a819fff95..6ccbe570e 100644 --- a/caiman/motion_correction.py +++ b/caiman/motion_correction.py @@ -167,9 +167,12 @@ def __init__(self, fname, min_mov=None, dview=None, max_shifts=(6, 6), niter_rig """ if 'ndarray' in str(type(fname)) or isinstance(fname, caiman.base.movies.movie): - logging.info('Creating file for motion correction "tmp_mov_mot_corr.hdf5"') - caiman.movie(fname).save('tmp_mov_mot_corr.hdf5') - fname = ['tmp_mov_mot_corr.hdf5'] + mc_tempfile = caiman.paths.fn_relocated('tmp_mov_mot_corr.hdf5') + if os.path.isfile(mc_tempfile): + os.remove(mc_tempfile) + logging.info(f"Creating file for motion correction: {mc_tempfile}") + caiman.movie(fname).save(mc_tempfile) + fname = [mc_tempfile] if not isinstance(fname, list): fname = [fname] @@ -2711,6 +2714,10 @@ def compute_metrics_motion_correction(fname, final_size_x, final_size_y, swap_di if play_flow and opencv: cv2.destroyAllWindows() + # FIXME: This generates a metrics dump potentially right next to the datafiles it was generated from; + # We will need to fix this in some future revision of the code, ideally returning the filename we used to the caller + # or abstracting the path handling logic into some kind of a policy-aware getpath function for specific uses. + # This should be fixed with future work to have separate runs have separate workdirs. np.savez(os.path.splitext(fname)[0] + '_metrics', flows=flows, norms=norms, correlations=correlations, smoothness=smoothness, tmpl=tmpl, smoothness_corr=smoothness_corr, img_corr=img_corr) return tmpl, correlations, flows, norms, smoothness @@ -3120,12 +3127,7 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 if base_name is None: base_name = os.path.splitext(os.path.split(fname)[1])[0] base_name = caiman.paths.fn_relocated(base_name) - fname_tot:Optional[str] = caiman.paths.memmap_frames_filename(base_name, dims, T, order) - if isinstance(fname, tuple): - fname_tot = os.path.join(os.path.split(fname[0])[0], fname_tot) - else: - fname_tot = os.path.join(os.path.split(fname)[0], fname_tot) np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=caiman.mmapping.prepare_shape(shape_mov), order=order) diff --git a/caiman/paths.py b/caiman/paths.py index ae5edbbe7..05abfb336 100644 --- a/caiman/paths.py +++ b/caiman/paths.py @@ -45,7 +45,7 @@ def get_tempdir() -> str: os.makedirs(temp_under_data) return temp_under_data -def fn_relocated(fn:str) -> str: +def fn_relocated(fn:str, force_temp:bool=False) -> str: """ If the provided filename does not contain any path elements, this returns what would be its absolute pathname as located in get_tempdir(). Otherwise it just returns what it is passed. @@ -53,10 +53,10 @@ def fn_relocated(fn:str) -> str: but if all they think about is filenames, they go under CaImAn's notion of its temporary dir. This is under the principle of "sensible defaults, but users can override them". """ - if not 'CAIMAN_NEW_TEMPFILE' in os.environ: # XXX We will ungate this in a future version of caiman - return fn - if str(os.path.basename(fn)) == str(fn): # No path stuff + if os.path.split(fn)[0] == '': # No path stuff return os.path.join(get_tempdir(), fn) + elif force_temp: + return os.path.join(get_tempdir(), os.path.split(fn)[1]) else: return fn @@ -124,4 +124,5 @@ def generate_fname_tot(base_name:str, dims:list[int], order:str) -> str: d1, d2, d3 = dims[0], dims[1], dims[2] ret = '_'.join([base_name, 'd1', str(d1), 'd2', str(d2), 'd3', str(d3), 'order', order]) ret = re.sub(r'(_)+', '_', ret) # Turn repeated underscores into just one - return ret + return fn_relocated(ret, force_temp=True) + diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 28f8718b5..b58e2e773 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -71,7 +71,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p ssub=2, tsub=2, p_ssub=1, p_tsub=1, method_init='greedy_roi', alpha_snmf=0.5, rf=None, stride=None, memory_fact=1, gnb=1, nb_patch=1, only_init_patch=False, method_deconvolution='oasis', n_pixels_per_process=4000, block_size_temp=5000, num_blocks_per_run_temp=20, - block_size_spat=5000, num_blocks_per_run_spat=20, + num_blocks_per_run_spat=20, check_nan=True, skip_refinement=False, normalize_init=True, options_local_NMF=None, minibatch_shape=100, minibatch_suff_stat=3, update_num_comps=True, rval_thr=0.9, thresh_fitness_delta=-20, @@ -86,7 +86,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p max_num_added=3, min_num_trial=2, thresh_CNN_noisy=0.5, fr=30, decay_time=0.4, min_SNR=2.5, ssub_B=2, init_iter=2, sniper_mode=False, use_peak_max=False, test_both=False, - expected_comps=500, max_merge_area=None, params=None): + expected_comps=500, params=None): """ Constructor of the CNMF method @@ -247,8 +247,6 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p init_iter: int, optional number of iterations for 1-photon imaging initialization - max_merge_area: int, optional - maximum area (in pixels) of merged components, used to determine whether to merge components during fitting process """ self.dview = dview @@ -272,7 +270,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p gnb=gnb, normalize_init=normalize_init, options_local_NMF=options_local_NMF, ring_size_factor=ring_size_factor, rolling_length=rolling_length, rolling_sum=rolling_sum, ssub=ssub, ssub_B=ssub_B, tsub=tsub, - block_size_spat=block_size_spat, num_blocks_per_run_spat=num_blocks_per_run_spat, + num_blocks_per_run_spat=num_blocks_per_run_spat, block_size_temp=block_size_temp, num_blocks_per_run_temp=num_blocks_per_run_temp, update_background_components=update_background_components, method_deconvolution=method_deconvolution, p=p, s_min=s_min, @@ -284,9 +282,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p n_refit=n_refit, num_times_comp_updated=num_times_comp_updated, simultaneously=simultaneously, sniper_mode=sniper_mode, test_both=test_both, thresh_CNN_noisy=thresh_CNN_noisy, thresh_fitness_delta=thresh_fitness_delta, thresh_fitness_raw=thresh_fitness_raw, thresh_overlap=thresh_overlap, - update_num_comps=update_num_comps, use_dense=use_dense, use_peak_max=use_peak_max, alpha_snmf=alpha_snmf, - max_merge_area=max_merge_area - ) + update_num_comps=update_num_comps, use_dense=use_dense, use_peak_max=use_peak_max, alpha_snmf=alpha_snmf) else: self.params = params params.set('patch', {'n_processes': n_processes}) @@ -479,7 +475,6 @@ def fit(self, images, indices=(slice(None), slice(None))): self.params.set('spatial', {'n_pixels_per_process': self.params.get('preprocess', 'n_pixels_per_process')}) logging.info('using ' + str(self.params.get('preprocess', 'n_pixels_per_process')) + ' pixels per process') - logging.info('using ' + str(self.params.get('spatial', 'block_size_spat')) + ' block_size_spat') logging.info('using ' + str(self.params.get('temporal', 'block_size_temp')) + ' block_size_temp') if self.params.get('patch', 'rf') is None: # no patches @@ -540,7 +535,7 @@ def fit(self, images, indices=(slice(None), slice(None))): logging.info('refinement...') if self.params.get('merging', 'do_merge'): logging.info('merging components ...') - self.merge_comps(Yr, mx=50, fast_merge=True, max_merge_area=self.params.get('merging', 'max_merge_area')) + self.merge_comps(Yr, mx=50, fast_merge=True) logging.info('Updating spatial ...') @@ -922,7 +917,7 @@ def update_spatial(self, Y, use_init=True, **kwargs): return self - def merge_comps(self, Y, mx=50, fast_merge=True, max_merge_area=None): + 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, \ @@ -933,8 +928,7 @@ def merge_comps(self, Y, mx=50, fast_merge=True, max_merge_area=None): self.params.get_group('spatial'), dview=self.dview, bl=self.estimates.bl, c1=self.estimates.c1, sn=self.estimates.neurons_sn, g=self.estimates.g, thr=self.params.get('merging', 'merge_thr'), mx=mx, - fast_merge=fast_merge, merge_parallel=self.params.get('merging', 'merge_parallel'), - max_merge_area=max_merge_area) + fast_merge=fast_merge, merge_parallel=self.params.get('merging', 'merge_parallel')) return self diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index c5889eb6e..1959a2e05 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -1235,7 +1235,7 @@ def deconvolve(self, params, dview=None, dff_flag=False): self.S_dff = np.stack([results[1][i] for i in order]) def merge_components(self, Y, params, mx=50, fast_merge=True, - dview=None, max_merge_area=None): + dview=None): """merges components """ # FIXME This method shares its name with a function elsewhere in the codebase (which it wraps) @@ -1247,8 +1247,7 @@ def merge_components(self, Y, params, mx=50, fast_merge=True, params.get_group('spatial'), dview=dview, bl=self.bl, c1=self.c1, sn=self.neurons_sn, g=self.g, thr=params.get('merging', 'merge_thr'), mx=mx, - fast_merge=fast_merge, merge_parallel=params.get('merging', 'merge_parallel'), - max_merge_area=max_merge_area) + fast_merge=fast_merge, merge_parallel=params.get('merging', 'merge_parallel')) def manual_merge(self, components, params): ''' merge a given list of components. The indices diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index 4a1528cd7..172f4e57b 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -288,7 +288,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter gSig = np.asarray(gSig, dtype=float) / ssub gSiz = np.round(np.asarray(gSiz) / ssub).astype(int) - if normalize_init is True: + if normalize_init: logging.info('Variance Normalization') if img is None: img = np.mean(Y, axis=-1) @@ -403,7 +403,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter Cin = np.empty((K, T), dtype=np.float32) center = [] - if normalize_init is True: + if normalize_init: if Ain.size > 0: Ain = Ain * np.reshape(img, (np.prod(d), -1), order='F') if sparse_b: @@ -657,7 +657,7 @@ def graphNMF(Y_ds, nr, max_iter_snmf=500, lambda_gnmf=1, C = mdl.fit_transform(yr).T A = mdl.components_.T W = caiman.source_extraction.cnmf.utilities.fast_graph_Laplacian_patches( - [np.reshape(m, [T, d], order='F').T, [], 'heat', SC_sigma, SC_thr, + [np.reshape(m, [T, d], order='F').T, [], SC_kernel, SC_sigma, SC_thr, SC_nnn, SC_normalize, SC_use_NN]) D = scipy.sparse.spdiags(W.sum(0), 0, W.shape[0], W.shape[0]) for it in range(max_iter_snmf): diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index baabb9723..5e6452232 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -19,7 +19,7 @@ def merge_components(Y, A, b, C, R, f, S, sn_pix, temporal_params, spatial_params, dview=None, thr=0.85, fast_merge=True, mx=1000, bl=None, c1=None, sn=None, g=None, - merge_parallel=False, max_merge_area=None) -> tuple[scipy.sparse.csc_matrix, np.ndarray, int, list, np.ndarray, float, float, float, float, list, np.ndarray]: + merge_parallel=False) -> tuple[scipy.sparse.csc_matrix, np.ndarray, int, list, np.ndarray, float, float, float, float, list, np.ndarray]: """ Merging of spatially overlapping components that have highly correlated temporal activity @@ -83,10 +83,6 @@ def merge_components(Y, A, b, C, R, f, S, sn_pix, temporal_params, merge_parallel: bool perform merging in parallel - max_merge_area: int - maximum area (in pixels) of merged components, - used to determine whether to merge - Returns: A: sparse matrix matrix of merged spatial components (d x K) diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index fefb1f100..2a1aede6c 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/caiman/source_extraction/cnmf/online_cnmf.py @@ -20,6 +20,7 @@ from math import sqrt from multiprocessing import cpu_count import numpy as np +import os from scipy.ndimage import percentile_filter from scipy.sparse import coo_matrix, csc_matrix, spdiags, hstack from scipy.stats import norm @@ -958,8 +959,9 @@ def initialize_online(self, model_LN=None, T=None): self.estimates = cnm.estimates else: - Y.save(caiman.paths.fn_relocated('init_file.hdf5')) - f_new = caiman.mmapping.save_memmap(['init_file.hdf5'], base_name='Yr', order='C', + temp_init_file = os.path.join(caiman.paths.get_tempdir(), 'init_file.hdf5') + Y.save(temp_init_file) + f_new = caiman.mmapping.save_memmap([temp_init_file], base_name='Yr', order='C', slices=[slice(0, opts['init_batch']), None, None]) Yrm, dims_, T_ = caiman.mmapping.load_memmap(f_new) diff --git a/caiman/source_extraction/cnmf/params.py b/caiman/source_extraction/cnmf/params.py index 595c8313d..7ca69fa04 100644 --- a/caiman/source_extraction/cnmf/params.py +++ b/caiman/source_extraction/cnmf/params.py @@ -29,7 +29,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), min_pnr=20, gnb=1, normalize_init=True, options_local_NMF=None, ring_size_factor=1.5, rolling_length=100, rolling_sum=True, ssub=2, ssub_B=2, tsub=2, - block_size_spat=5000, num_blocks_per_run_spat=20, + num_blocks_per_run_spat=20, block_size_temp=5000, num_blocks_per_run_temp=20, update_background_components=True, method_deconvolution='oasis', p=2, s_min=None, @@ -97,12 +97,6 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), last_commit: str hash of last commit in the caiman repo. Pleaes do not override this. - mmap_F: list[str] - paths to F-order memory mapped files after motion correction - - mmap_C: str - path to C-order memory mapped file after motion correction - CNMFParams.patch (these control how the data is divided into patches): border_pix: int, default: 0 Number of pixels to exclude around each border. @@ -142,7 +136,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), If list, it should be a list of two elements corresponding to the height and width of patches skip_refinement: bool, default: False - Whether to skip refinement of components (deprecated?) + Whether to skip refinement of components p_ssub: float, default: 2 Spatial downsampling factor @@ -157,7 +151,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), check_nan: bool, default: True whether to check for NaNs - compute_g': bool, default: False + compute_g: bool, default: False whether to estimate global time constant include_noise: bool, default: False @@ -191,7 +185,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), K: int, default: 30 number of components to be found (per patch or whole FOV depending on whether rf=None) - SC_kernel: {'heat', 'cos', binary'}, default: 'heat' + SC_kernel: {'heat', 'cos', 'binary'}, default: 'heat' kernel for graph affinity matrix SC_sigma: float, default: 1 @@ -288,9 +282,6 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), temporal downsampling factor CNMFParams.spatial (these control how the algorithms handle spatial components): - block_size_spat : int, default: 5000 - Number of pixels to process at the same time for dot product. Reduce if you face memory problems - dist: float, default: 3 expansion factor of ellipse @@ -318,7 +309,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), number of pixels to be processed by each worker nb: int, default: 1 - number of global background components + number of global background components. Do not set this directly; modify it in init. normalize_yyt_one: bool, default: True Whether to normalize the C and A matrices so that diag(C*C.T) = 1 during update spatial @@ -361,15 +352,12 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), optimize_g: bool, default: False flag for optimizing time constants - memory_efficient: - (undocumented) - method_deconvolution: 'oasis'|'cvxpy'|'oasis', default: 'oasis' method for solving the constrained deconvolution problem ('oasis','cvx' or 'cvxpy') if method cvxpy, primary and secondary (if problem unfeasible for approx solution) nb: int, default: 1 - number of global background components + number of global background components. Do not set this directly; modify it in init. noise_method: 'mean'|'median'|'logmexp', default: 'mean' PSD averaging method for computing the noise std @@ -402,9 +390,6 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), merge_parallel: bool, default: False Perform merging in parallel - max_merge_area: int or None, default: None - maximum area (in pixels) of merged components, used to determine whether to merge components during fitting process - CNMFParams.quality (these control how quality of traces are evaluated): SNR_lowest: float, default: 0.5 minimum required trace SNR. Traces with SNR below this will get rejected @@ -677,9 +662,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'dxy': dxy, 'var_name_hdf5': var_name_hdf5, 'caiman_version': pkg_resources.get_distribution('caiman').version, - 'last_commit': None, - 'mmap_F': None, - 'mmap_C': None + 'last_commit': None } self.patch = { @@ -716,7 +699,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), } self.init = { - 'K': k, # number of components, + 'K': k, # number of components, 'SC_kernel': 'heat', # kernel for graph affinity matrix 'SC_sigma' : 1, # std for SC kernel 'SC_thr': 0, # threshold for affinity matrix @@ -756,7 +739,6 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), } self.spatial = { - 'block_size_spat': block_size_spat, # number of pixels to parallelize residual computation ** DECREASE IF MEMORY ISSUES 'dist': 3, # expansion factor of ellipse 'expandCore': iterate_structure(generate_binary_structure(2, 1), 2).astype(int), # Flag to extract connected components (might want to turn to False for dendritic imaging) @@ -793,7 +775,6 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), # number of autocovariance lags to be considered for time constant estimation 'lags': 5, 'optimize_g': False, # flag for optimizing time constants - 'memory_efficient': False, # method for solving the constrained deconvolution problem ('oasis','cvx' or 'cvxpy') # if method cvxpy, primary and secondary (if problem unfeasible for approx # solution) solvers to be used with cvxpy, can be 'ECOS','SCS' or 'CVXOPT' @@ -811,8 +792,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), self.merging = { 'do_merge': do_merge, 'merge_thr': merge_thresh, - 'merge_parallel': False, - 'max_merge_area': max_merge_area + 'merge_parallel': False } self.quality = { @@ -824,7 +804,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'rval_lowest': -1, # minimum accepted space correlation 'rval_thr': rval_thr, # space correlation threshold 'use_cnn': True, # use CNN based classifier - 'use_ecc': False, # flag for eccentricity based filtering + 'use_ecc': False, # flag for eccentricity based filtering (2D only) 'max_ecc': 3 } @@ -957,8 +937,6 @@ def check_consistency(self): if self.patch['rf'] is not None: if np.any(np.array(self.patch['rf']) <= self.init['gSiz'][0]): logging.warning(f"Changing rf from {self.patch['rf']} to {2 * self.init['gSiz'][0]} because the constraint rf > gSiz was not satisfied.") -# if self.motion['gSig_filt'] is None: -# self.motion['gSig_filt'] = self.init['gSig'] if self.init['nb'] <= 0 and (self.patch['nb_patch'] != self.init['nb'] or self.patch['low_rank_background'] is not None): logging.warning(f"gnb={self.init['nb']}, hence setting keys nb_patch and low_rank_background in group patch automatically.") @@ -1145,7 +1123,7 @@ def change_params(self, params_dict, allow_legacy:bool=True, warn_unused:bool=Tr consumed = {} # Keep track of what parameters in params_dict were used to set something in params (just for legacy API) nagged_once = False # So we don't nag people multiple times in the same call for paramkey in params_dict: - if paramkey in list(self.__dict__.keys()): # Proper pathed part + if paramkey in list(self.__dict__.keys()) and isinstance(params_dict[paramkey], dict): # Handle proper pathed part. Latter half of the conditional is because of scoped keys with the same name as categories, because we apparently have those. ring_CNN is an example. cat_handle = getattr(self, paramkey) for k, v in params_dict[paramkey].items(): if k == 'nb' and paramkey != 'init': diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 38451ea87..f8b1a833d 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -40,7 +40,7 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, se=np.ones((3, 3), dtype=int), ss=np.ones((3, 3), dtype=int), nb=1, method_ls='lasso_lars', update_background_components=True, - low_rank_background=True, block_size_spat=1000, + low_rank_background=True, num_blocks_per_run_spat=20): """update spatial footprints and background through Basis Pursuit Denoising diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 549a4d01b..b5bdde791 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -114,9 +114,6 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N You should start the cluster (install ipyparallel and then type ipcluster -n 6, where 6 is the number of processes). - memory_efficient: Bool - whether or not to optimize for memory usage (longer running times). necessary with very large datasets - kwargs: dict all parameters passed to constrained_foopsi except bl,c1,g,sn (see documentation). Some useful parameters are @@ -287,9 +284,6 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, You should start the cluster (install ipyparallel and then type ipcluster -n 6, where 6 is the number of processes). - memory_efficient: Bool - whether or not to optimize for memory usage (longer running times). necessary with very large datasets - **kwargs: dict all parameters passed to constrained_foopsi except bl,c1,g,sn (see documentation). Some useful parameters are diff --git a/caiman/tests/comparison_humans.py b/caiman/tests/comparison_humans.py index ea12f1fd4..5b955e26b 100644 --- a/caiman/tests/comparison_humans.py +++ b/caiman/tests/comparison_humans.py @@ -347,7 +347,6 @@ 'p': global_params['p'], }, 'spatial': { - 'block_size_spat': block_size, 'nb': global_params['gnb'], 'num_blocks_per_run_spat': num_blocks_per_run, 'n_pixels_per_process': n_pixels_per_process, diff --git a/caiman/tests/test_motion_correction.py b/caiman/tests/test_motion_correction.py index 98a72d5dd..af56a3050 100644 --- a/caiman/tests/test_motion_correction.py +++ b/caiman/tests/test_motion_correction.py @@ -126,8 +126,7 @@ def gen_data(D=2, noise=.01, T=300, framerate=30, firerate=2., motion=True): def _test_motion_correct_rigid(D): Y, C, S, A, centers, dims, shifts = gen_data(D) - fname = 'testMovie.tif' - cm.movie(Y).save(fname) + fname = cm.movie(Y).save('testMovie.tif') params_dict = { 'motion': { 'border_nan': True, diff --git a/caiman/tests/test_sbx.py b/caiman/tests/test_sbx.py index c54459222..d8609edaa 100644 --- a/caiman/tests/test_sbx.py +++ b/caiman/tests/test_sbx.py @@ -4,6 +4,7 @@ import numpy.testing as npt import os import tifffile +import tracemalloc import caiman as cm from caiman.paths import caiman_datadir @@ -90,37 +91,61 @@ def test_load_subind(): npt.assert_array_equal(data_3d_plane0_3d[:, :, :, 0], data_3d_plane0_2d) -def test_sbx_to_tif(): - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest1.tif') +def test_load_efficiency(): + # Make sure that when loading, excess copies are not being made and + # data outside subindices are not being loaded into memory file_2d = os.path.join(TESTDATA_PATH, '2d_sbx.sbx') - data_2d_from_sbx = cm.load(file_2d) - sbx_utils.sbx_to_tif(file_2d, fileout=tif_file) - data_2d_from_tif = cm.load(tif_file) - npt.assert_array_almost_equal(data_2d_from_sbx, data_2d_from_tif, - err_msg='Data do not match when loaded from .sbx vs. .tif') + tracemalloc.start() + data_2d_sliced = sbx_utils.sbxread(file_2d, subindices=(slice(None), slice(None, None, 2))) + curr_mem, peak_mem = tracemalloc.get_traced_memory() + assert peak_mem / curr_mem < 1.1, 'Too much memory allocated when loading' + del data_2d_sliced + tracemalloc.stop() - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest2.tif') - file_3d = os.path.join(TESTDATA_PATH, '3d_sbx_1.sbx') - data_3d_from_sbx = cm.load(file_3d, is3D=True) - sbx_utils.sbx_to_tif(file_3d, fileout=tif_file) - data_3d_from_tif = cm.load(tif_file, is3D=True) - npt.assert_array_almost_equal(data_3d_from_sbx, data_3d_from_tif, - err_msg='3D data do not match when loaded from .sbx vs. .tif') - # make sure individual planes are not saved as 3D (i.e. RGB) - Y = tifffile.TiffFile(tif_file).series[0] - assert Y.shape == SHAPE_3D, 'Shape of data in tif file is wrong' - if Y[0].shape != SHAPE_3D[2:]: - if Y[0].shape == SHAPE_3D[1:]: - assert False, 'Tif "plane" is 3-dimensional (i.e., has channel dimension)' - else: - assert False, 'Shape of tif plane is wrong' - - # with subindices - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest3.tif') - subinds = (slice(0, None, 2), [0, 1, 3], slice(None)) - sbx_utils.sbx_to_tif(file_2d, fileout=tif_file, subindices=subinds) - sub_data_from_tif = cm.load(tif_file) - npt.assert_array_almost_equal(data_2d_from_sbx[subinds_to_ix(subinds, data_2d_from_sbx.shape)], sub_data_from_tif) + +def test_sbx_to_tif(): + tif_filename = os.path.join(caiman_datadir(), 'temp', 'from_sbx.tif') + tif_file = None + + try: + file_2d = os.path.join(TESTDATA_PATH, '2d_sbx.sbx') + data_2d_from_sbx = cm.load(file_2d) + sbx_utils.sbx_to_tif(file_2d, fileout=tif_filename) + data_2d_from_tif = cm.load(tif_filename) + npt.assert_array_almost_equal(data_2d_from_sbx, data_2d_from_tif, + err_msg='Data do not match when loaded from .sbx vs. .tif') + + file_3d = os.path.join(TESTDATA_PATH, '3d_sbx_1.sbx') + data_3d_from_sbx = cm.load(file_3d, is3D=True) + sbx_utils.sbx_to_tif(file_3d, fileout=tif_filename) + data_3d_from_tif = cm.load(tif_filename, is3D=True) + npt.assert_array_almost_equal(data_3d_from_sbx, data_3d_from_tif, + err_msg='3D data do not match when loaded from .sbx vs. .tif') + # make sure individual planes are not saved as 3D (i.e. RGB) + with tifffile.TiffFile(tif_filename) as tif_file: + Y = tif_file.series[0] + assert Y.shape == SHAPE_3D, 'Shape of data in tif file is wrong' + if Y[0].shape != SHAPE_3D[2:]: + if Y[0].shape == SHAPE_3D[1:]: + assert False, 'Tif "plane" is 3-dimensional (i.e., has channel dimension)' + else: + assert False, 'Shape of tif plane is wrong' + + # with subindices + subinds = (slice(0, None, 2), [0, 1, 3], slice(None)) + sbx_utils.sbx_to_tif(file_2d, fileout=tif_filename, subindices=subinds) + sub_data_from_tif = cm.load(tif_filename) + npt.assert_array_almost_equal(data_2d_from_sbx[subinds_to_ix(subinds, data_2d_from_sbx.shape)], sub_data_from_tif) + + # with plane + sbx_utils.sbx_to_tif(file_3d, fileout=tif_filename, plane=0) + plane_data_from_tif = cm.load(tif_filename) + npt.assert_array_almost_equal(data_3d_from_sbx[:, :, :, 0], plane_data_from_tif) + + finally: + # cleanup + if os.path.isfile(tif_filename): + os.remove(tif_filename) # with plane tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest4.tif') @@ -129,37 +154,41 @@ def test_sbx_to_tif(): npt.assert_array_almost_equal(data_3d_from_sbx[:, :, :, 0], plane_data_from_tif) def test_sbx_chain_to_tif(): - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest5.tif') - file_3d_1 = os.path.join(TESTDATA_PATH, '3d_sbx_1.sbx') - data_3d_1 = sbx_utils.sbxread(file_3d_1) - file_3d_2 = os.path.join(TESTDATA_PATH, '3d_sbx_2.sbx') - data_3d_2 = sbx_utils.sbxread(file_3d_2) - - # normal chain - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest6.tif') - sbx_utils.sbx_chain_to_tif([file_3d_1, file_3d_2], fileout=tif_file) - data_chain_tif = cm.load(tif_file, is3D=True) - data_chain_gt = np.concatenate([data_3d_1, data_3d_2], axis=0) - npt.assert_array_almost_equal(data_chain_tif, data_chain_gt, - err_msg='Tif from chain does not match expected') - - # matching subindices - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest7.tif') - sbx_utils.sbx_chain_to_tif([file_3d_1, file_3d_2], fileout=tif_file, - subindices=(slice(None), slice(0, None, 2))) - data_chain_tif = cm.load(tif_file, is3D=True) - data_chain_gt = data_chain_gt[:, ::2] - npt.assert_array_almost_equal(data_chain_tif, data_chain_gt, - err_msg='Tif from chain with subindices does not match expected') + tif_filename = os.path.join(caiman_datadir(), 'temp', 'from_sbx.tif') + try: + file_3d_1 = os.path.join(TESTDATA_PATH, '3d_sbx_1.sbx') + data_3d_1 = sbx_utils.sbxread(file_3d_1) + file_3d_2 = os.path.join(TESTDATA_PATH, '3d_sbx_2.sbx') + data_3d_2 = sbx_utils.sbxread(file_3d_2) + + # normal chain + sbx_utils.sbx_chain_to_tif([file_3d_1, file_3d_2], fileout=tif_filename) + data_chain_tif = cm.load(tif_filename, is3D=True) + data_chain_gt = np.concatenate([data_3d_1, data_3d_2], axis=0) + npt.assert_array_almost_equal(data_chain_tif, data_chain_gt, + err_msg='Tif from chain does not match expected') + + # matching subindices + sbx_utils.sbx_chain_to_tif([file_3d_1, file_3d_2], fileout=tif_filename, + subindices=(slice(None), slice(0, None, 2))) + data_chain_tif = cm.load(tif_filename, is3D=True) + data_chain_gt = data_chain_gt[:, ::2] + npt.assert_array_almost_equal(data_chain_tif, data_chain_gt, + err_msg='Tif from chain with subindices does not match expected') + + # non-matching subindices with compatible shapes + subinds_1 = (slice(None), [0, 1, 3], slice(0, None, 2), [0, 2]) + subinds_2 = (slice(1, None), [-4, -2, -1], slice(1, None, 2), [1, 3]) + sbx_utils.sbx_chain_to_tif([file_3d_1, file_3d_2], fileout=tif_filename, + subindices=[subinds_1, subinds_2]) + data_chain_tif = cm.load(tif_filename, is3D=True) + data_chain_gt = np.concatenate([data_3d_1[subinds_to_ix(subinds_1, data_3d_1.shape)], + data_3d_2[subinds_to_ix(subinds_2, data_3d_2.shape)]], axis=0) + npt.assert_array_almost_equal(data_chain_tif, data_chain_gt, + err_msg='Tif from chain with non-matching subindices does not match expected') + + finally: + # cleanup + if os.path.isfile(tif_filename): + os.remove(tif_filename) - # non-matching subindices with compatible shapes - tif_file = os.path.join(caiman_datadir(), 'temp', 'sbxtest8.tif') - subinds_1 = (slice(None), [0, 1, 3], slice(0, None, 2), [0, 2]) - subinds_2 = (slice(1, None), [-4, -2, -1], slice(1, None, 2), [1, 3]) - sbx_utils.sbx_chain_to_tif([file_3d_1, file_3d_2], fileout=tif_file, - subindices=[subinds_1, subinds_2]) - data_chain_tif = cm.load(tif_file, is3D=True) - data_chain_gt = np.concatenate([data_3d_1[subinds_to_ix(subinds_1, data_3d_1.shape)], - data_3d_2[subinds_to_ix(subinds_2, data_3d_2.shape)]], axis=0) - npt.assert_array_almost_equal(data_chain_tif, data_chain_gt, - err_msg='Tif from chain with non-matching subindices does not match expected') diff --git a/caiman/tests/test_toydata.py b/caiman/tests/test_toydata.py index f55316f83..2a0f2a8e0 100644 --- a/caiman/tests/test_toydata.py +++ b/caiman/tests/test_toydata.py @@ -44,7 +44,6 @@ def pipeline(D): gSig=[2, 2, 2][:D], p=1, n_pixels_per_process=np.prod(dims), - block_size_spat=np.prod(dims), block_size_temp=np.prod(dims)) params.spatial['thr_method'] = 'nrg' params.spatial['extract_cc'] = False diff --git a/caiman/utils/sbx_utils.py b/caiman/utils/sbx_utils.py index b8ff03254..f349b9d0b 100644 --- a/caiman/utils/sbx_utils.py +++ b/caiman/utils/sbx_utils.py @@ -80,6 +80,7 @@ def sbxread(filename: str, subindices: Optional[FileSubindices] = slice(None), c def sbx_to_tif(filename: str, fileout: Optional[str] = None, subindices: Optional[FileSubindices] = slice(None), + bigtiff: Optional[bool] = True, imagej: bool = False, to32: bool = False, channel: Optional[int] = None, plane: Optional[int] = None, chunk_size: int = 1000): """ Convert a single .sbx file to .tif format @@ -95,6 +96,9 @@ def sbx_to_tif(filename: str, fileout: Optional[str] = None, subindices: Optiona which frames to read (defaults to all) if a tuple of non-scalars, specifies slices of up to 4 dimensions in the order (frame, Y, X, Z). + to32: bool + whether to save in float32 format (default is to keep as uint16) + channel: int | None which channel to save (required if data has >1 channel) @@ -106,31 +110,17 @@ def sbx_to_tif(filename: str, fileout: Optional[str] = None, subindices: Optiona how many frames to load into memory at once (None = load the whole thing) """ # Check filenames - basename, ext = os.path.splitext(filename) - if ext == '.sbx': - filename = basename - if fileout is None: + basename, ext = os.path.splitext(filename) + if ext == '.sbx': + filename = basename fileout = filename + '.tif' - else: - # Add a '.tif' extension if not already present - extension = os.path.splitext(fileout)[1].lower() - if extension not in ['.tif', '.tiff', '.btf']: - fileout = fileout + '.tif' if subindices is None: subindices = slice(None) - # Find the shape of the file we need - save_shape = _get_output_shape(filename, subindices)[0] - if plane is not None: - if len(save_shape) < 4: - raise Exception('Plane cannot be specified for 2D data') - save_shape = save_shape[:3] - - # Open tif as a memmap and copy to it - memmap_tif = tifffile.memmap(fileout, shape=save_shape, dtype='uint16', photometric='MINISBLACK') - _sbxread_helper(filename, subindices=subindices, channel=channel, out=memmap_tif, plane=plane, chunk_size=chunk_size) + sbx_chain_to_tif([filename], fileout, [subindices], bigtiff=bigtiff, imagej=imagej, to32=to32, + channel=channel, plane=plane, chunk_size=chunk_size) def sbx_chain_to_tif(filenames: list[str], fileout: str, subindices: Optional[ChainSubindices] = slice(None), @@ -149,11 +139,8 @@ def sbx_chain_to_tif(filenames: list[str], fileout: str, subindices: Optional[Ch see subindices for sbx_to_tif can specify separate subindices for each file if nested 2 levels deep; X, Y, and Z sizes must match for all files after indexing. - - to32: bool - whether to save in float32 format (default is to keep as uint16) - channel, plane, chunk_size: see sbx_to_tif + to32, channel, plane, chunk_size: see sbx_to_tif """ if subindices is None: subindices = slice(None) @@ -187,33 +174,37 @@ def sbx_chain_to_tif(filenames: list[str], fileout: str, subindices: Optional[Ch channel = 0 elif any(shape[0] <= channel for shape in all_shapes): raise Exception('Not all files have the requested channel') - + # Allocate empty tif file with the final shape (do this first to ensure any existing file is overwritten) common_shape = tuple(map(int, all_shapes_out[0, 1:])) - Ns = list(map(int, all_shapes_out[:, 0])) - N = sum(Ns) - save_shape = (N,) + common_shape + all_n_frames_out = list(map(int, all_shapes_out[:, 0])) + n_frames_out = sum(all_n_frames_out) + save_shape = (n_frames_out,) + common_shape if plane is not None: if len(save_shape) < 4: raise Exception('Plane cannot be specified for 2D data') save_shape = save_shape[:3] + # Add a '.tif' extension if not already present extension = os.path.splitext(fileout)[1].lower() if extension not in ['.tif', '.tiff', '.btf']: fileout = fileout + '.tif' dtype = np.float32 if to32 else np.uint16 + # Make the file first so we can pass in bigtiff and imagej options; otherwise could create using tifffile.memmap directly tifffile.imwrite(fileout, data=None, shape=save_shape, bigtiff=bigtiff, imagej=imagej, - dtype=dtype, photometric='MINISBLACK') + dtype=dtype, photometric='MINISBLACK', align=tifffile.TIFF.ALLOCATIONGRANULARITY) # Now convert each file tif_memmap = tifffile.memmap(fileout, series=0) offset = 0 - for filename, subind, file_N in zip(filenames, subindices, Ns): + for filename, subind, file_N in zip(filenames, subindices, all_n_frames_out): _sbxread_helper(filename, subindices=subind, channel=channel, out=tif_memmap[offset:offset+file_N], plane=plane, chunk_size=chunk_size) offset += file_N + del tif_memmap # important to make sure file is closed (on Windows) + def sbx_shape(filename: str, info: Optional[dict] = None) -> tuple[int, int, int, int, int]: """ @@ -282,7 +273,7 @@ def sbx_shape(filename: str, info: Optional[dict] = None) -> tuple[int, int, int else: info['max_idx'] = filesize / info['bytesPerBuffer'] * factor - 1 - N = info['max_idx'] + 1 # Last frame + n_frames = info['max_idx'] + 1 # Last frame # Determine whether we are looking at a z-stack # Only consider optotune z-stacks - knobby schedules have too many possibilities and @@ -291,9 +282,9 @@ def sbx_shape(filename: str, info: Optional[dict] = None) -> tuple[int, int, int n_planes = info['otparam'][2] else: n_planes = 1 - N //= n_planes + n_frames //= n_planes - x = (int(info['nChan']), int(info['sz'][1]), int(info['recordsPerBuffer']), int(n_planes), int(N)) + x = (int(info['nChan']), int(info['sz'][1]), int(info['recordsPerBuffer']), int(n_planes), int(n_frames)) return x @@ -440,69 +431,54 @@ def _sbxread_helper(filename: str, subindices: FileSubindices = slice(None), cha raise Exception('Invalid scanbox metadata') save_shape, subindices = _get_output_shape(data_shape, subindices) - N = save_shape[0] + n_frames_out = save_shape[0] if plane is not None: if len(save_shape) < 4: raise Exception('Plane cannot be specified for 2D data') save_shape = save_shape[:3] - frame_stride = _interpret_subindices(subindices[0], n_frames)[1] - if out is not None and out.shape != save_shape: raise Exception('Existing memmap is the wrong shape to hold loaded data') - # Open File - with open(filename + '.sbx') as sbx_file: - def get_chunk(n: int): - # Note: SBX files store the values strangely, it's necessary to invert each uint16 value to get the correct ones - chunk = np.invert(np.fromfile(sbx_file, dtype='uint16', count=frame_size//2 * n)) - chunk = np.reshape(chunk, data_shape[:4] + (n,), order='F') # (chans, X, Y, Z, frames) - chunk = np.transpose(chunk, (0, 4, 2, 1, 3))[channel] # to (frames, Y, X, Z) - if not is3D: # squeeze out singleton plane dim - chunk = np.squeeze(chunk, axis=3) - chunk = chunk[(slice(None),) + np.ix_(*subindices[1:])] - if is3D and len(save_shape) == 3: # squeeze out plane dim after selecting - chunk = chunk[:, :, :, plane] - return chunk - - if frame_stride == 1: - sbx_file.seek(subindices[0][0] * frame_size, 0) - - if chunk_size is None: - # load a contiguous block all at once - if out is None: - # avoid allocating an extra array just to copy into it - out = get_chunk(N) - else: - out[:] = get_chunk(N) - else: - # load a contiguous block of data, but chunk it to not use too much memory at once - if out is None: - out = np.empty(save_shape, dtype=np.uint16) - - n_remaining = N - offset = 0 - while n_remaining > 0: - this_chunk_size = min(n_remaining, chunk_size) - chunk = get_chunk(this_chunk_size) - out[offset:offset+this_chunk_size] = chunk - n_remaining -= this_chunk_size - offset += this_chunk_size - + # Read from .sbx file, using memmap to avoid loading until necessary + sbx_mmap = np.memmap(filename + '.sbx', mode='r', dtype='uint16', shape=data_shape, order='F') + sbx_mmap = np.transpose(sbx_mmap, (0, 4, 2, 1, 3)) # to (chans, frames, Y, X, Z) + sbx_mmap = sbx_mmap[channel] + if not is3D: # squeeze out singleton plane dim + sbx_mmap = sbx_mmap[..., 0] + elif plane is not None: # select plane relative to subindices + sbx_mmap = sbx_mmap[..., subindices[-1][plane]] + subindices = subindices[:-1] + inds = np.ix_(*subindices) + + if chunk_size is None: + # load a contiguous block all at once + chunk_size = n_frames_out + elif out is None: + # Pre-allocate destination when loading in chunks + out = np.empty(save_shape, dtype=np.uint16) + + n_remaining = n_frames_out + offset = 0 + while n_remaining > 0: + this_chunk_size = min(n_remaining, chunk_size) + # Note: important to copy the data here instead of making a view, + # so the memmap can be closed (achieved by advanced indexing) + chunk = sbx_mmap[(inds[0][offset:offset+this_chunk_size],) + inds[1:]] + # Note: SBX files store the values strangely, it's necessary to invert each uint16 value to get the correct ones + np.invert(chunk, out=chunk) # avoid copying, may be large + + if out is None: + out = chunk # avoid copying when loading all data else: - # load one frame at a time, sorting indices for fastest access - if out is None: - out = np.empty(save_shape, dtype=np.uint16) - - for counter, (k, ind) in enumerate(sorted(zip(subindices[0], range(N)))): - if counter % 100 == 0: - logging.debug(f'Reading iteration: {k}') - sbx_file.seek(k * frame_size, 0) - out[ind] = get_chunk(1) - - if isinstance(out, np.memmap): - out.flush() - + out[offset:offset+this_chunk_size] = chunk + n_remaining -= this_chunk_size + offset += this_chunk_size + + del sbx_mmap # Important to close file (on Windows) + + if isinstance(out, np.memmap): + out.flush() return out @@ -520,7 +496,10 @@ def _interpret_subindices(subindices: DimSubindices, dim_extent: int) -> tuple[I f'(requested up to {subindices.stop})') else: iterable_elements = subindices - skip = 0 + if isinstance(subindices, range): + skip = subindices.step + else: + skip = 0 return iterable_elements, skip diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index d1a5832e4..39e0db78f 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -528,7 +528,7 @@ def nb_view_patches3d(Y_r, A, C, dims, image_type='mean', Yr=None, if max_projection: if image_type == 'corr': - tmp = [(local_correlations( + tmp = [(caiman.summary_images.local_correlations( Yr[index_permut].reshape(dims + (-1,), order='F'))[:, ::-1]).max(i) for i in range(3)] diff --git a/demos/README.txt b/demos/README.txt new file mode 100644 index 000000000..4229e2837 --- /dev/null +++ b/demos/README.txt @@ -0,0 +1,27 @@ +The demos serve a few purposes: +1) They demonstrate use of the Caiman software with some sample data either + bundled or automatically fetchable +2) They act as a template for you to write your own analysis +3) They can be directly used if they already do almost exactly what you need + to do with Caiman. + +There are two types of demos, each in their own folder: +* Jupyter demos are meant to be run interactively, in Jupyter or potentially + one of its variants or something else compatible with Jupyter notebooks. + They run in a cell-oriented way, allowing you to step through each part of + the Caiman pipeline, seeing results and possibly changing things (values of + parameters, files to work on, similar) between them. +* CLI demos are meant to be run from the commandline of your operating system + (Powershell or CMD on Windows, Bash/Zsh inside a Terminal on Linux or OSX). + When running Caiman this way, you might or might not want to see + visualisations, and the arguments you provide to these demos control their + behaviour. + +The CLI demos used to have a different intended use case, also being meant for +cell-based interactive use, but in an IPython-based IDE like Spyder; this usage +is no longer supported (use something compatible with Jupyter notebook formats +for that). The earlier form of the demos is available in our use_cases repo in +this directory: + +https://github.com/flatironinstitute/caiman_use_cases/tree/main/use_cases/old_cli_demos + diff --git a/demos/general/demo_OnACID.py b/demos/general/demo_OnACID.py index 32bd44b11..1951e6375 100755 --- a/demos/general/demo_OnACID.py +++ b/demos/general/demo_OnACID.py @@ -1,92 +1,66 @@ #!/usr/bin/env python -# -*- coding: utf-8 -*- """ Basic demo for the CaImAn Online algorithm (OnACID) using CNMF initialization. It demonstrates the construction of the params and online_cnmf objects and the fit function that is used to run the algorithm. For a more complete demo check the script demo_OnACID_mesoscope.py - -@author: jfriedrich & epnev """ -#%% -from IPython import get_ipython + +import argparse +#import code import logging +import matplotlib import numpy as np import os +try: + cv2.setNumThreads(0) +except: + pass + import caiman as cm from caiman.source_extraction import cnmf as cnmf from caiman.paths import caiman_datadir -try: - if __IPYTHON__: - print("Detected iPython") - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: - pass -#%% -# Set up the logger; change this if you like. -# You can log to a file using the filename parameter, or make the output more or less -# verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR - -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]"\ - "[%(process)d] %(message)s", - level=logging.WARNING) - # filename="/tmp/caiman.log" -# %% def main(): - pass # For compatibility between running under an IDE and the CLI - - # %% load data - fname = [os.path.join(caiman_datadir(), 'example_movies', 'demoMovie.tif')] - - # %% set up some parameters - fr = 10 # frame rate (Hz) - decay_time = .75 # approximate length of transient event in seconds - gSig = [6, 6] # expected half size of neurons - p = 1 # order of AR indicator dynamics - min_SNR = 1 # minimum SNR for accepting candidate components - thresh_CNN_noisy = 0.65 # CNN threshold for candidate components - gnb = 2 # number of background components - init_method = 'cnmf' # initialization method - - # set up CNMF initialization parameters - init_batch = 400 # number of frames for initialization - patch_size = 32 # size of patch - stride = 3 # amount of overlap between patches - K = 4 # max number of components in each patch - params_dict = {'fr': fr, - 'fnames': fname, - 'decay_time': decay_time, - 'gSig': gSig, - 'p': p, - 'min_SNR': min_SNR, - 'nb': gnb, - 'init_batch': init_batch, - 'init_method': init_method, - 'rf': patch_size//2, - 'stride': stride, - 'sniper_mode': True, - 'thresh_CNN_noisy': thresh_CNN_noisy, - 'K': K} - opts = cnmf.params.CNMFParams(params_dict=params_dict) - - # %% fit with online object + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + opts = cnmf.params.CNMFParams(params_from_file=cfg.configfile) + + if cfg.input is not None: + opts.change_params({"data": {"fnames": cfg.input}}) + + if not opts.data['fnames']: # Set neither by CLI arg nor through JSON, so use default data + fnames = [os.path.join(caiman_datadir(), 'example_movies', 'demoMovie.tif')] + opts.change_params({"data": {"fnames": fnames}}) + + # If you want to break into an interactive console session, move and uncomment this wherever you want in the code + # (and uncomment the code import at the top) + #code.interact(local=dict(globals(), **locals()) ) + + # fit with online object cnm = cnmf.online_cnmf.OnACID(params=opts) cnm.fit_online() - # %% plot contours - logging.info('Number of components:' + str(cnm.estimates.A.shape[-1])) - Cn = cm.load(fname[0], subindices=slice(0,500)).local_correlations(swap_dim=False) + # plot contours + logging.info(f"Number of components: {cnm.estimates.A.shape[-1]}") + Cn = cm.load(fnames[0], subindices=slice(0,500)).local_correlations(swap_dim=False) cnm.estimates.plot_contours(img=Cn) - # %% pass through the CNN classifier with a low threshold (keeps clearer neuron shapes and excludes processes) + # pass through the CNN classifier with a low threshold (keeps clearer neuron shapes and excludes processes) use_CNN = True if use_CNN: # threshold for CNN classifier @@ -94,11 +68,18 @@ def main(): cnm.estimates.evaluate_components_CNN(opts) cnm.estimates.plot_contours(img=Cn, idx=cnm.estimates.idx_components) - # %% plot results + # plot results cnm.estimates.view_components(img=Cn, idx=cnm.estimates.idx_components) + if not cfg.no_play: + matplotlib.pyplot.show(block=True) + +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate basic Caiman Online functionality with CNMF initialization") + parser.add_argument("--configfile", default=os.path.join(caiman_datadir(), 'demos', 'general', 'params_demo_OnACID.json'), help="JSON Configfile for Caiman parameters") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() -# %% -# This is to mask the differences between running this demo in IDE -# versus from the CLI -if __name__ == "__main__": - main() +######## +main() diff --git a/demos/general/demo_OnACID_mesoscope.py b/demos/general/demo_OnACID_mesoscope.py index cf9bde48f..620a8b1b9 100755 --- a/demos/general/demo_OnACID_mesoscope.py +++ b/demos/general/demo_OnACID_mesoscope.py @@ -1,5 +1,4 @@ #!/usr/bin/env python -# -*- coding: utf-8 -*- """ Complete pipeline for online processing using CaImAn Online (OnACID). @@ -11,27 +10,21 @@ can be passed as options. A plot of the processing time for the various steps of the algorithm is also included. -@author: Eftychios Pnevmatikakis @epnev - -Special thanks to Andreas Tolias and his lab at Baylor College of Medicine +Thanks to Andreas Tolias and his lab at Baylor College of Medicine for sharing the data used in this demo. """ +import argparse import glob -from IPython import get_ipython +import logging +import matplotlib import numpy as np import os -import logging import matplotlib.pyplot as plt try: - if __IPYTHON__: - print("Detected iPython") - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: + cv2.setNumThreads(0) +except: pass import caiman as cm @@ -40,92 +33,56 @@ from caiman.utils.utils import download_demo -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]"\ - "[%(process)d] %(message)s", - level=logging.INFO) - -# %% def main(): - pass # For compatibility between running under IDE and the CLI - - # %% download and list all files to be processed - - # folder inside ./example_movies where files will be saved - fld_name = 'Mesoscope' - fnames = [] - fnames.append(download_demo('Tolias_mesoscope_1.hdf5', fld_name)) - fnames.append(download_demo('Tolias_mesoscope_2.hdf5', fld_name)) - fnames.append(download_demo('Tolias_mesoscope_3.hdf5', fld_name)) - - # your list of files should look something like this - logging.info(fnames) - - # %% Set up some parameters - - fr = 15 # frame rate (Hz) - decay_time = 0.5 # approximate length of transient event in seconds - gSig = (3, 3) # expected half size of neurons - p = 1 # order of AR indicator dynamics - min_SNR = 1 # minimum SNR for accepting new components - ds_factor = 1 # spatial downsampling factor (increases speed but may lose some fine structure) - gnb = 2 # number of background components - gSig = tuple(np.ceil(np.array(gSig) / ds_factor).astype('int')) # recompute gSig if downsampling is involved - mot_corr = True # flag for online motion correction - pw_rigid = False # flag for pw-rigid motion correction (slower but potentially more accurate) - max_shifts_online = np.ceil(10.).astype('int') # maximum allowed shift during motion correction - sniper_mode = True # use a CNN to detect new neurons (o/w space correlation) - rval_thr = 0.9 # soace correlation threshold for candidate components - # set up some additional supporting parameters needed for the algorithm - # (these are default values but can change depending on dataset properties) - init_batch = 200 # number of frames for initialization (presumably from the first file) - K = 2 # initial number of components - epochs = 1 # number of passes over the data - show_movie = False # show the movie as the data gets processed - - params_dict = {'fnames': fnames, - 'fr': fr, - 'decay_time': decay_time, - 'gSig': gSig, - 'p': p, - 'min_SNR': min_SNR, - 'rval_thr': rval_thr, - 'ds_factor': ds_factor, - 'nb': gnb, - 'motion_correct': mot_corr, - 'init_batch': init_batch, - 'init_method': 'bare', - 'normalize': True, - 'sniper_mode': sniper_mode, - 'K': K, - 'epochs': epochs, - 'max_shifts_online': max_shifts_online, - 'pw_rigid': pw_rigid, - 'dist_shape_update': True, - 'min_num_trial': 10, - 'show_movie': show_movie} - opts = cnmf.params.CNMFParams(params_dict=params_dict) - - # %% fit online + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + opts = cnmf.params.CNMFParams(params_from_file=cfg.configfile) + + if cfg.input is not None and len(cfg.input) != 0: # CLI arg can override all other settings for fnames, although other data-centric commands still must match source data + opts.change_params({'data': {'fnames': cfg.input}}) + + if not opts.data['fnames']: # Set neither by CLI arg nor through JSON, so use default data + fld_name = 'Mesoscope' + print("Downloading demos...") + fnames = [download_demo('Tolias_mesoscope_1.hdf5', fld_name), + download_demo('Tolias_mesoscope_2.hdf5', fld_name), + download_demo('Tolias_mesoscope_3.hdf5', fld_name)] + opts.change_params({'data': {'fnames': fnames}}) + + opts.change_params({'online': {'show_movie': cfg.show_movie}}) # Override from CLI + + print(f"Params: {opts.to_json()}") + # fit online cnm = cnmf.online_cnmf.OnACID(params=opts) cnm.fit_online() - # %% plot contours (this may take time) - logging.info('Number of components: ' + str(cnm.estimates.A.shape[-1])) + # plot contours (this may take time) + logging.info(f"Number of components: {cnm.estimates.A.shape[-1]}") images = cm.load(fnames) Cn = images.local_correlations(swap_dim=False, frames_per_chunk=500) cnm.estimates.plot_contours(img=Cn, display_numbers=False) - # %% view components + # view components cnm.estimates.view_components(img=Cn) - # %% plot timing performance (if a movie is generated during processing, timing + # plot timing performance (if a movie is generated during processing, timing # will be severely over-estimated) T_motion = 1e3*np.array(cnm.t_motion) T_detect = 1e3*np.array(cnm.t_detect) T_shapes = 1e3*np.array(cnm.t_shapes) - T_track = 1e3*np.array(cnm.t_online) - T_motion - T_detect - T_shapes + T_track = 1e3*np.array(cnm.t_online) - T_motion - T_detect - T_shapes plt.figure() plt.stackplot(np.arange(len(T_motion)), T_motion, T_track, T_detect, T_shapes) plt.legend(labels=['motion', 'tracking', 'detect', 'shapes'], loc=2) @@ -133,10 +90,8 @@ def main(): plt.xlabel('Frame #') plt.ylabel('Processing time [ms]') - #%% RUN IF YOU WANT TO VISUALIZE THE RESULTS (might take time) - c, dview, n_processes = \ - cm.cluster.setup_cluster(backend='multiprocessing', n_processes=None, - single_thread=False) + # Prepare result visualisations (might take time) + c, dview, n_processes = cm.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) if opts.online['motion_correct']: shifts = cnm.estimates.shifts[-cnm.estimates.C.shape[-1]:] if not opts.motion['pw_rigid']: @@ -157,26 +112,26 @@ def main(): Yr, dims, T = cm.load_memmap(memmap_file) images = np.reshape(Yr.T, [T] + list(dims), order='F') - min_SNR = 2 # peak SNR for accepted components (if above this, accept) - rval_thr = 0.85 # space correlation threshold (if above this, accept) - use_cnn = True # use the CNN classifier - min_cnn_thr = 0.99 # if cnn classifier predicts below this value, reject - cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected - - cnm.params.set('quality', {'min_SNR': min_SNR, - 'rval_thr': rval_thr, - 'use_cnn': use_cnn, - 'min_cnn_thr': min_cnn_thr, - 'cnn_lowest': cnn_lowest}) + # evaluate_components uses parameters from the 'quality' category cnm.estimates.evaluate_components(images, cnm.params, dview=dview) cnm.estimates.Cn = Cn - cnm.save(os.path.splitext(fnames[0])[0]+'_results.hdf5') + cnm.save(os.path.splitext(fnames[0])[0] + '_results.hdf5') dview.terminate() - -#%% -# This is to mask the differences between running this demo in IDE -# versus from the CLI -if __name__ == "__main__": - main() \ No newline at end of file + if not cfg.no_play: + matplotlib.pyplot.show(block=True) + +def handle_args(): + parser = argparse.ArgumentParser(description="Full OnACID Caiman demo") + parser.add_argument("--configfile", default=os.path.join(caiman_datadir(), 'demos', 'general', 'params_demo_OnACID_mesoscope.json'), help="JSON Configfile for Caiman parameters") + parser.add_argument("--show_movie", action="store_true", help="Display results as movie") + parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") + parser.add_argument("--cluster_nproc", type=int, default=None, help="Override automatic selection of number of workers to use") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() + +######## +main() \ No newline at end of file diff --git a/demos/general/demo_behavior.py b/demos/general/demo_behavior.py index 9547d4e39..19c46224e 100755 --- a/demos/general/demo_behavior.py +++ b/demos/general/demo_behavior.py @@ -1,95 +1,89 @@ #!/usr/bin/env python -# -*- coding: utf-8 -*- """ -Created on Sun Aug 6 11:43:39 2017 +Demonstrate Caiman functions relating to behavioral experiments and optical flow -@author: agiovann +This demo requires a GUI; it does not make sense to run it noninteractively. """ +# If this demo crashes right after making it past the initial play of the mouse movie, +# you may need to switch your build of opencv+libopencv+py-opencv to a qt6 build of the +# same (in conda, these have a qt6_ prefix to their build id). + +import argparse import cv2 -from IPython import get_ipython import logging import numpy as np -import matplotlib.pyplot as plt +import matplotlib.pyplot as pl try: cv2.setNumThreads(0) except: pass -try: - if __IPYTHON__: - print("Detected iPython") - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: - pass - import caiman as cm from caiman.behavior import behavior from caiman.utils.utils import download_demo -#%% -# Set up the logger; change this if you like. -# You can log to a file using the filename parameter, or make the output more or less -# verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR - -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s", - # filename="/tmp/caiman.log", - level=logging.WARNING) -#%% def main(): - pass # For compatibility between running under IDE and the CLI + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + if cfg.input is None: + # If no input is specified, use sample data, downloading if necessary + fnames = [download_demo('demo_behavior.h5')] + else: + fnames = cfg.input + # If you prefer to hardcode filenames, you could do something like this: + # fnames = ["/path/to/myfile1.avi", "/path/to/myfile2.avi"] + + pl.ion() - #%% - plt.ion() - - fname = [u'demo_behavior.h5'] - if fname[0] in ['demo_behavior.h5']: - # TODO: todocument - fname = [download_demo(fname[0])] # TODO: todocument - m = cm._load_behavior(fname[0]) + m = cm._load_behavior(fnames[0]) - #%% load, rotate and eliminate useless pixels - m = m.transpose([0, 2, 1]) - m = m[:, 150:, :] + # load, rotate and eliminate useless pixels + m = m.transpose([0, 2, 1]) # XXX Does this really work outside this dataset? + m = m[:, 150:, :] # TODO adopt some syntax for clipping, or make this optional and tell the user to clip before running - #%% visualize movie + # visualize movie m.play() - #%% select interesting portion of the FOV (draw a polygon on the figure that pops up, when done press enter) - # TODO: Put the message below into the image + # select interesting portion of the FOV (draw a polygon on the figure that pops up, when done press enter) print("Please draw a polygon delimiting the ROI on the image that will be displayed after the image; press enter when done") mask = np.array(behavior.select_roi(np.median(m[::100], 0), 1)[0], np.float32) - #%% n_components = 4 # number of movement looked for resize_fact = 0.5 # for computational efficiency movies are downsampled # number of standard deviations above mean for the magnitude that are considered enough to measure the angle in polar coordinates num_std_mag_for_angle = .6 only_magnitude = False # if onlu interested in factorizing over the magnitude - method_factorization = 'dict_learn' # could also use nmf + method_factorization = 'nmf' # number of iterations for the dictionary learning algorithm (Marial et al, 2010) max_iter_DL = -30 spatial_filter_, time_trace_, of_or = cm.behavior.behavior.extract_motor_components_OF(m, n_components, mask=mask, - resize_fact=resize_fact, only_magnitude=only_magnitude, verbose=True, method_factorization='dict_learn', max_iter_DL=max_iter_DL) + resize_fact=resize_fact, only_magnitude=only_magnitude, verbose=True, method_factorization='nmf', max_iter_DL=max_iter_DL) + - #%% mags, dircts, dircts_thresh, spatial_masks_thrs = cm.behavior.behavior.extract_magnitude_and_angle_from_OF( spatial_filter_, time_trace_, of_or, num_std_mag_for_angle=num_std_mag_for_angle, sav_filter_size=3, only_magnitude=only_magnitude) - - #%% + idd = 0 - axlin = plt.subplot(n_components, 2, 2) + axlin = pl.subplot(n_components, 2, 2) for mag, dirct, spatial_filter in zip(mags, dircts_thresh, spatial_filter_): - plt.subplot(n_components, 2, 1 + idd * 2) + pl.subplot(n_components, 2, 1 + idd * 2) min_x, min_y = np.min(np.where(mask), 1) spfl = spatial_filter @@ -99,19 +93,23 @@ def main(): mask[min_x:max_x, min_y:max_y] = spfl mask[mask < np.nanpercentile(spfl, 70)] = np.nan - plt.imshow(m[0], cmap='gray') - plt.imshow(mask, alpha=.5) - plt.axis('off') + pl.imshow(m[0], cmap='gray') + pl.imshow(mask, alpha=.5) + pl.axis('off') - axelin = plt.subplot(n_components, 2, 2 + idd * 2, sharex=axlin) - plt.plot(mag / 10, 'k') + axelin = pl.subplot(n_components, 2, 2 + idd * 2, sharex=axlin) + pl.plot(mag / 10, 'k') dirct[mag < 0.5 * np.std(mag)] = np.nan - plt.plot(dirct, 'r-', linewidth=2) + pl.plot(dirct, 'r-', linewidth=2) idd += 1 + pl.show(block=True) + +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate behavioural/optic flow functions") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() -#%% -# This is to mask the differences between running this demo in an IDE -# versus from the CLI -if __name__ == "__main__": - main() \ No newline at end of file +######## +main() diff --git a/demos/general/demo_caiman_basic.py b/demos/general/demo_caiman_basic.py index b5ce1e7e0..255e5fb00 100755 --- a/demos/general/demo_caiman_basic.py +++ b/demos/general/demo_caiman_basic.py @@ -1,9 +1,8 @@ #!/usr/bin/env python -# -*- coding: utf-8 -*- """ -Basic stripped-down demo for running the CNMF source extraction algorithm with -CaImAn and evaluation the components. The analysis can be run either in the +Basic demo for running the CNMF source extraction algorithm with +Caiman and evaluation the components. The analysis can be run either in the whole FOV or in patches. For a complete pipeline (including motion correction) check demo_pipeline.py @@ -12,13 +11,11 @@ This demo is designed to be run in an IDE or from a CLI; its plotting functions are tailored for that environment. -@authors: @agiovann and @epnev - """ +import argparse import cv2 import glob -from IPython import get_ipython import logging import numpy as np import os @@ -28,144 +25,110 @@ except: pass -try: - if __IPYTHON__: - print("Detected iPython") - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') -except NameError: - pass - - -import caiman as cm +import caiman from caiman.paths import caiman_datadir from caiman.source_extraction.cnmf import cnmf as cnmf from caiman.source_extraction.cnmf import params as params from caiman.summary_images import local_correlations_movie_offline +from caiman.utils.utils import download_demo -# %% -# Set up the logger; change this if you like. -# You can log to a file using the filename parameter, or make the output more or less -# verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]"\ - "[%(process)d] %(message)s", - level=logging.WARNING) - # filename="/tmp/caiman.log", - -#%% def main(): - pass # For compatibility between running under an IDE or a CLI - - # %% start a cluster - - c, dview, n_processes =\ - cm.cluster.setup_cluster(backend='multiprocessing', n_processes=None, - single_thread=False) - - # %% set up some parameters - fnames = [os.path.join(caiman_datadir(), 'example_movies', 'demoMovie.tif')] - # file(s) to be analyzed - is_patches = True # flag for processing in patches or not - fr = 10 # approximate frame rate of data - decay_time = 5.0 # length of transient - - if is_patches: # PROCESS IN PATCHES AND THEN COMBINE - rf = 10 # half size of each patch - stride = 4 # overlap between patches - K = 4 # number of components in each patch - else: # PROCESS THE WHOLE FOV AT ONCE - rf = None # setting these parameters to None - stride = None # will run CNMF on the whole FOV - K = 30 # number of neurons expected (in the whole FOV) - - gSig = [6, 6] # expected half size of neurons - merge_thresh = 0.80 # merging threshold, max correlation allowed - p = 2 # order of the autoregressive system - gnb = 2 # global background order - - params_dict = {'fnames': fnames, - 'fr': fr, - 'decay_time': decay_time, - 'rf': rf, - 'stride': stride, - 'K': K, - 'gSig': gSig, - 'merge_thr': merge_thresh, - 'p': p, - 'nb': gnb} - - opts = params.CNMFParams(params_dict=params_dict) - - # %% Now RUN CaImAn Batch (CNMF) + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + opts = params.CNMFParams(params_from_file=cfg.configfile) + + if cfg.input is not None: # CLI arg can override all other settings for fnames, although other data-centric commands still must match source data + opts.change_params({'data': {'fnames': cfg.input}}) + + if not opts.data['fnames']: # Set neither by CLI arg nor through JSON, so use default data + # If no input is specified, use sample data, downloading if necessary + fnames = [os.path.join(caiman_datadir(), 'example_movies', 'demoMovie.tif')] # file(s) to be analyzed + opts.change_params({'data': {'fnames': fnames}}) + + # start a cluster for parallel processing + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) + + + # Run CaImAn Batch (CNMF) cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) cnm = cnm.fit_file() - # %% plot contour plots of components + # plot contour plots of components Cns = local_correlations_movie_offline(fnames[0], remove_baseline=True, swap_dim=False, window=1000, stride=1000, winSize_baseline=100, quantil_min_baseline=10, dview=dview) Cn = Cns.max(axis=0) - cnm.estimates.plot_contours(img=Cn) + Cn[np.isnan(Cn)] = 0 # Convert NaNs to zero + if not cfg.no_play: + cnm.estimates.plot_contours(img=Cn) - # %% load memory mapped file - Yr, dims, T = cm.load_memmap(cnm.mmap_file) + # load memory mapped file + Yr, dims, T = caiman.load_memmap(cnm.mmap_file) images = np.reshape(Yr.T, [T] + list(dims), order='F') - # %% refit + # refit cnm2 = cnm.refit(images, dview=dview) - # %% COMPONENT EVALUATION - # the components are evaluated in three ways: + # Component Evaluation + # Components are evaluated in three ways: # a) the shape of each component must be correlated with the data # b) a minimum peak SNR is required over the length of a transient # c) each shape passes a CNN based classifier (this will pick up only neurons # and filter out active processes) - min_SNR = 2 # peak SNR for accepted components (if above this, accept) - rval_thr = 0.85 # space correlation threshold (if above this, accept) - use_cnn = True # use the CNN classifier - min_cnn_thr = 0.99 # if cnn classifier predicts below this value, reject - cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected - - cnm2.params.set('quality', {'min_SNR': min_SNR, - 'rval_thr': rval_thr, - 'use_cnn': use_cnn, - 'min_cnn_thr': min_cnn_thr, - 'cnn_lowest': cnn_lowest}) - - cnm2.estimates.evaluate_components(images, cnm2.params, dview=dview) - # %% visualize selected and rejected components - cnm2.estimates.plot_contours(img=Cn, idx=cnm2.estimates.idx_components) + cnm2.estimates.evaluate_components(images, opts, dview=dview) - # %% visualize selected components - cnm2.estimates.view_components(images, idx=cnm2.estimates.idx_components, img=Cn) + if not cfg.no_play: + # visualize selected and rejected components + cnm2.estimates.plot_contours(img=Cn, idx=cnm2.estimates.idx_components) + # visualize selected components + cnm2.estimates.view_components(images, idx=cnm2.estimates.idx_components, img=Cn) - #%% only select high quality components (destructive) + # only select high quality components (destructive) # cnm2.estimates.select_components(use_object=True) # cnm2.estimates.plot_contours(img=Cn) - #%% save results + # save results cnm2.estimates.Cn = Cn - cnm2.save(cnm2.mmap_file[:-4]+'hdf5') - - # %% play movie with results (original, reconstructed, amplified residual) - cnm2.estimates.play_movie(images, magnification=4); - - # %% STOP CLUSTER and clean up log files - cm.stop_server(dview=dview) - - log_files = glob.glob('Yr*_LOG_*') - for log_file in log_files: - os.remove(log_file) - -# %% -# This is to mask the differences between running this demo in IDE -# versus from the CLI -if __name__ == "__main__": - main() \ No newline at end of file + cnm2.save(cnm2.mmap_file[:-4] + 'hdf5') # FIXME use the path library to do this the right way + + # play movie with results (original, reconstructed, amplified residual) + if not cfg.no_play: + cnm2.estimates.play_movie(images, magnification=4); + + # Stop the cluster and clean up log files + caiman.stop_server(dview=dview) + + if not cfg.keep_logs: + log_files = glob.glob('Yr*_LOG_*') + for log_file in log_files: + os.remove(log_file) + +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate basic Caiman functionality") + parser.add_argument("--configfile", default=os.path.join(caiman_datadir(), 'demos', 'general', 'params_demo_caiman_basic_with_patches.json'), help="JSON Configfile for Caiman parameters") + parser.add_argument("--keep_logs", action="store_true", help="Keep temporary logfiles") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") + parser.add_argument("--cluster_nproc", type=int, default=None, help="Override automatic selection of number of workers to use") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() + +######## +main() diff --git a/demos/general/demo_pipeline.py b/demos/general/demo_pipeline.py index 4d50ecc94..028b9c55a 100755 --- a/demos/general/demo_pipeline.py +++ b/demos/general/demo_pipeline.py @@ -2,7 +2,7 @@ """ Complete demo pipeline for processing two photon calcium imaging data using the -CaImAn batch algorithm. The processing pipeline included motion correction, +Caiman batch algorithm. The processing pipeline included motion correction, source extraction and deconvolution. The demo shows how to construct the params, MotionCorrect and cnmf objects and call the relevant functions. You can also run a large part of the pipeline with a single method (cnmf.fit_file) @@ -13,16 +13,13 @@ This demo pertains to two photon data. For a complete analysis pipeline for one photon microendoscopic data see demo_pipeline_cnmfE.py - -copyright GNU General Public License v2.0 -authors: @agiovann and @epnev """ -#%% + +import argparse import cv2 import glob -from IPython import get_ipython import logging -import matplotlib.pyplot as plt +import matplotlib import numpy as np import os @@ -31,164 +28,88 @@ except: pass -try: - if __IPYTHON__: - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: - pass - - -import caiman as cm +import caiman from caiman.motion_correction import MotionCorrect from caiman.source_extraction.cnmf import cnmf as cnmf from caiman.source_extraction.cnmf import params as params -from caiman.utils.utils import download_demo from caiman.summary_images import local_correlations_movie_offline +from caiman.utils.utils import download_demo + -# %% -# Set up the logger (optional); change this if you like. -# You can log to a file using the filename parameter, or make the output more -# or less verbose by setting level to logging.DEBUG, logging.INFO, -# logging.WARNING, or logging.ERROR -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]"\ - "[%(process)d] %(message)s", - level=logging.WARNING) - -#%% def main(): - pass # For compatibility between running in an IDE and the CLI - - #%% Select file(s) to be processed (download if not present) - fnames = ['Sue_2x_3000_40_-46.tif'] # filename to be processed - if fnames[0] in ['Sue_2x_3000_40_-46.tif', 'demoMovie.tif']: - fnames = [download_demo(fnames[0])] - - #%% First setup some parameters for data and motion correction - - # dataset dependent parameters - fr = 30 # imaging rate in frames per second - decay_time = 0.4 # length of a typical transient in seconds - dxy = (2., 2.) # spatial resolution in x and y in (um per pixel) - # note the lower than usual spatial resolution here - max_shift_um = (12., 12.) # maximum shift in um - patch_motion_um = (100., 100.) # patch size for non-rigid correction in um - - # motion correction parameters - pw_rigid = True # flag to select rigid vs pw_rigid motion correction - # maximum allowed rigid shift in pixels - max_shifts = [int(a/b) for a, b in zip(max_shift_um, dxy)] - # start a new patch for pw-rigid motion correction every x pixels - strides = tuple([int(a/b) for a, b in zip(patch_motion_um, dxy)]) - # overlap between patches (size of patch in pixels: strides+overlaps) - overlaps = (24, 24) - # maximum deviation allowed for patch with respect to rigid shifts - max_deviation_rigid = 3 - - mc_dict = { - 'fnames': fnames, - 'fr': fr, - 'decay_time': decay_time, - 'dxy': dxy, - 'pw_rigid': pw_rigid, - 'max_shifts': max_shifts, - 'strides': strides, - 'overlaps': overlaps, - 'max_deviation_rigid': max_deviation_rigid, - 'border_nan': 'copy' - } - - opts = params.CNMFParams(params_dict=mc_dict) - - # %% play the movie (optional) + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + # First set up some parameters for data and motion correction + opts = params.CNMFParams(params_from_file=cfg.configfile) + + if cfg.input is not None: + opts.change_params({"data": {"fnames": cfg.input}}) + if not opts.data['fnames']: # Set neither by CLI arg nor through JSON, so use default data + fnames = [download_demo('Sue_2x_3000_40_-46.tif')] + opts.change_params({"data": {"fnames": fnames}}) + + m_orig = caiman.load_movie_chain(opts.data['fnames']) + + # play the movie (optional) # playing the movie using opencv. It requires loading the movie in memory. # To close the video press q - display_images = False - if display_images: - m_orig = cm.load_movie_chain(fnames) + if not cfg.no_play: ds_ratio = 0.2 moviehandle = m_orig.resize(1, 1, ds_ratio) moviehandle.play(q_max=99.5, fr=60, magnification=2) - # %% start a cluster for parallel processing - c, dview, n_processes = cm.cluster.setup_cluster( - backend='multiprocessing', n_processes=None, single_thread=False) + # start a cluster for parallel processing + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) - # %%% MOTION CORRECTION + # Motion Correction # first we create a motion correction object with the specified parameters - mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion')) + mc = MotionCorrect(opts.data['fnames'], dview=dview, **opts.get_group('motion')) # note that the file is not loaded in memory - # %% Run (piecewise-rigid motion) correction using NoRMCorre + # Run (piecewise-rigid motion) correction using NoRMCorre mc.motion_correct(save_movie=True) - # %% compare with original movie - if display_images: - m_orig = cm.load_movie_chain(fnames) - m_els = cm.load(mc.mmap_file) + # compare with original movie + if not cfg.no_play: + m_orig = caiman.load_movie_chain(opts.data['fnames']) + m_els = caiman.load(mc.mmap_file) ds_ratio = 0.2 - moviehandle = cm.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie, + moviehandle = caiman.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie, m_els.resize(1, 1, ds_ratio)], axis=2) moviehandle.play(fr=60, q_max=99.5, magnification=2) # press q to exit - # %% MEMORY MAPPING + # MEMORY MAPPING border_to_0 = 0 if mc.border_nan == 'copy' else mc.border_to_0 # you can include the boundaries of the FOV if you used the 'copy' option # during motion correction, although be careful about the components near # the boundaries # memory map the file in order 'C' - fname_new = cm.save_memmap(mc.mmap_file, base_name='memmap_', order='C', + fname_new = caiman.save_memmap(mc.mmap_file, base_name='memmap_', order='C', border_to_0=border_to_0) # exclude borders # now load the file - Yr, dims, T = cm.load_memmap(fname_new) + Yr, dims, T = caiman.load_memmap(fname_new) images = np.reshape(Yr.T, [T] + list(dims), order='F') # load frames in python format (T x X x Y) - # %% restart cluster to clean up memory - cm.stop_server(dview=dview) - c, dview, n_processes = cm.cluster.setup_cluster( - backend='multiprocessing', n_processes=None, single_thread=False) - - # %% parameters for source extraction and deconvolution - p = 1 # order of the autoregressive system - gnb = 2 # number of global background components - merge_thr = 0.85 # merging threshold, max correlation allowed - rf = 15 - # half-size of the patches in pixels. e.g., if rf=25, patches are 50x50 - stride_cnmf = 6 # amount of overlap between the patches in pixels - K = 4 # number of components per patch - gSig = [4, 4] # expected half size of neurons in pixels - # initialization method (if analyzing dendritic data using 'sparse_nmf') - method_init = 'greedy_roi' - ssub = 2 # spatial subsampling during initialization - tsub = 2 # temporal subsampling during initialization - - # parameters for component evaluation - opts_dict = {'fnames': fnames, - 'p': p, - 'fr': fr, - 'nb': gnb, - 'rf': rf, - 'K': K, - 'gSig': gSig, - 'stride': stride_cnmf, - 'method_init': method_init, - 'rolling_sum': True, - 'merge_thr': merge_thr, - 'n_processes': n_processes, - 'only_init': True, - 'ssub': ssub, - 'tsub': tsub} - - opts.change_params(params_dict=opts_dict); - - # %% RUN CNMF ON PATCHES + # restart cluster to clean up memory + caiman.stop_server(dview=dview) + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) + + # RUN CNMF ON PATCHES # First extract spatial and temporal components on patches and combine them # for this step deconvolution is turned off (p=0). If you want to have # deconvolution within each patch change params.patch['p_patch'] to a @@ -198,86 +119,82 @@ def main(): cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) cnm = cnm.fit(images) - # %% Alternate way to run above pipeline using a single method (optional) - # you can also perform the motion correction plus cnmf fitting steps - # simultaneously after defining your parameters object using - # cnm1 = cnmf.CNMF(n_processes, params=opts, dview=dview) - # cnm1.fit_file(motion_correct=True) - - # %% plot contours of found components + # plot contours of found components Cns = local_correlations_movie_offline(mc.mmap_file[0], - remove_baseline=True, window=1000, stride=1000, + remove_baseline=True, + window=1000, stride=1000, winSize_baseline=100, quantil_min_baseline=10, dview=dview) Cn = Cns.max(axis=0) Cn[np.isnan(Cn)] = 0 - cnm.estimates.plot_contours(img=Cn) - plt.title('Contour plots of found components'); + if not cfg.no_play: + cnm.estimates.plot_contours(img=Cn) - #%% save results + # save results cnm.estimates.Cn = Cn - cnm.save(fname_new[:-5]+'_init.hdf5') + cnm.save(fname_new[:-5] + '_init.hdf5') # FIXME - # %% RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution + # RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution cnm2 = cnm.refit(images, dview=dview) - # %% COMPONENT EVALUATION - # the components are evaluated in three ways: + # Component Evaluation + # Components are evaluated in three ways: # a) the shape of each component must be correlated with the data # b) a minimum peak SNR is required over the length of a transient # c) each shape passes a CNN based classifier - min_SNR = 2 # signal to noise ratio for accepting a component - rval_thr = 0.85 # space correlation threshold for accepting a component - cnn_thr = 0.99 # threshold for CNN based classifier - cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected - - cnm2.params.set('quality', {'decay_time': decay_time, - 'min_SNR': min_SNR, - 'rval_thr': rval_thr, - 'use_cnn': True, - 'min_cnn_thr': cnn_thr, - 'cnn_lowest': cnn_lowest}) + cnm2.estimates.evaluate_components(images, cnm2.params, dview=dview) - # %% PLOT COMPONENTS - cnm2.estimates.plot_contours(img=Cn, idx=cnm2.estimates.idx_components) + if not cfg.no_play: + # Plot Components + cnm2.estimates.plot_contours(img=Cn, idx=cnm2.estimates.idx_components) - # %% VIEW TRACES (accepted and rejected) - if display_images: + # View Traces (accepted and rejected) cnm2.estimates.view_components(images, img=Cn, idx=cnm2.estimates.idx_components) cnm2.estimates.view_components(images, img=Cn, idx=cnm2.estimates.idx_components_bad) - - #%% update object with selected components (optional) + + # update object with selected components (optional) cnm2.estimates.select_components(use_object=True) - #%% Extract DF/F values + # Extract DF/F values cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250) - #%% Show final traces - cnm2.estimates.view_components(img=Cn) + # Show final traces + if not cfg.no_play: + cnm2.estimates.view_components(img=Cn) - #%% cnm2.estimates.Cn = Cn cnm2.save(cnm2.mmap_file[:-4] + 'hdf5') - #%% reconstruct denoised movie (press q to exit) - if display_images: + # reconstruct denoised movie (press q to exit) + if not cfg.no_play: cnm2.estimates.play_movie(images, q_max=99.9, gain_res=2, magnification=2, bpx=border_to_0, include_bck=False) # background not shown - #%% STOP CLUSTER and clean up log files - cm.stop_server(dview=dview) - log_files = glob.glob('*_LOG_*') - for log_file in log_files: - os.remove(log_file) - -# %% -# This is to mask the differences between running this demo in an IDE -# versus from the CLI -if __name__ == "__main__": - main() + # Stop the cluster and clean up log files + caiman.stop_server(dview=dview) + matplotlib.pyplot.show(block=True) + + if not cfg.keep_logs: + log_files = glob.glob('*_LOG_*') + for log_file in log_files: + os.remove(log_file) + +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate 2P Pipeline using batch algorithm") + parser.add_argument("--configfile", default=os.path.join(caiman.paths.caiman_datadir(), 'demos', 'general', 'params_demo_pipeline.json'), help="JSON Configfile for Caiman parameters") + parser.add_argument("--keep_logs", action="store_true", help="Keep temporary logfiles") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") + parser.add_argument("--cluster_nproc", type=int, default=None, help="Override automatic selection of number of workers to use") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() + +######## +main() diff --git a/demos/general/demo_pipeline_NWB.py b/demos/general/demo_pipeline_NWB.py old mode 100644 new mode 100755 index f1335adcd..956f632b1 --- a/demos/general/demo_pipeline_NWB.py +++ b/demos/general/demo_pipeline_NWB.py @@ -4,12 +4,15 @@ This script follows closely the demo_pipeline.py script but uses the Neurodata Without Borders (NWB) file format for loading the input and saving the output. It is meant as an example on how to use NWB files with CaImAn. -authors: @agiovann and @epnev + +If you provide a filename in the config json, it must be an nwb file. """ -#%% + +import argparse import cv2 +from datetime import datetime +from dateutil.tz import tzlocal import glob -from IPython import get_ipython import logging import matplotlib.pyplot as plt import numpy as np @@ -20,271 +23,227 @@ except: pass -try: - if __IPYTHON__: - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: - pass - - -from datetime import datetime -from dateutil.tz import tzlocal - -import caiman as cm +import caiman from caiman.motion_correction import MotionCorrect +from caiman.paths import caiman_datadir, get_tempdir from caiman.source_extraction.cnmf import cnmf as cnmf from caiman.source_extraction.cnmf import params as params from caiman.utils.utils import download_demo -from caiman.paths import caiman_datadir - -# %% -# Set up the logger (optional); change this if you like. -# You can log to a file using the filename parameter, or make the output more -# or less verbose by setting level to logging.DEBUG, logging.INFO, -# logging.WARNING, or logging.ERROR -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]"\ - "[%(process)d] %(message)s", - level=logging.WARNING) - -#%% + +# TODO: Various code here isn't careful enough with paths and may produce files in +# places alongside their sources, rather than in temporary directories. Come back +# later and fix this once we've worked out where everything should go. + def main(): - pass # For compatibility between running under an IDE and the CLI + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + opts = params.CNMFParams(params_from_file=cfg.configfile) + if cfg.input is not None: + opts.change_params({"data": {"fnames": cfg.input}}) + # TODO throw error if it is not an nwb file + if not opts.data['fnames']: # Set neither by CLI arg nor through JSON, so use default data + # This demo will either accept a nwb file or it'll convert an existing movie into one + # default path, no input was provide, so we'll take the "Sue" sample data, convert it to a nwb, and then ingest it. + # The convert target is the original fn, with the extension swapped to nwb. + fr = float(15) # imaging rate in frames per second + pre_convert = download_demo('Sue_2x_3000_40_-46.tif') + convert_target = os.path.join(get_tempdir(), 'Sue_2x_3000_40_-46.nwb') + + orig_movie = caiman.load(pre_convert, fr=fr) + # save file in NWB format with various additional info + orig_movie.save(convert_target, + sess_desc="test", + identifier="demo 1", + imaging_plane_description='single plane', + emission_lambda=520.0, indicator='GCAMP6f', + location='parietal cortex', + experimenter='Sue Ann Koay', lab_name='Tank Lab', + institution='Princeton U', + experiment_description='Experiment Description', + session_id='Session 1', + var_name_hdf5='TwoPhotonSeries') + fnames = [convert_target] + opts.change_params({'data': {'fnames': fnames, 'var_name_hdf5': 'TwoPhotonSeries'}}) + elif cfg.input[0].endswith('.nwb'): + if os.path.isfile(cfg.input[0]): + # We were handed at least one nwb file and it exists either right here or as an absolute path + for fn in cfg.input: + if not os.path.isfile(fn): + raise Exception(f"Could not find input file {fn}") + if not fn.endswith('.nwb'): + raise Exception(f"Cannot mix nwb and other file formats in input with this demo") + fnames = cfg.input + elif os.path.isfile(os.path.join(caiman_datadir(), cfg.input[0])): + # We were handed at least one nwb file and it exists, but in the caiman_datadir() + # So repeat the logic above except look in caiman_datadir() + for fn in cfg.input: + if not os.path.isfile(os.path.join(caiman_datadir(), fn)): + raise Exception(f"Could not find input file {fn}") + if not fn.endswith('.nwb'): + raise Exception(f"Cannot mix nwb and other file formats in input with this demo") + fnames = list(map(lambda n: os.path.join(caiman_datadir, n), cfg.input)) + + else: # Someone mentioned an nwb file but we can't find it! + raise Exception(f"Could not find referenced input file {cfg.input[0]}") + else: + raise Exception("If you're providing your own data to this demo, you must provide it in nwb format or pre-convert it yourself") - #%% Select file(s) to be processed (download if not present) - fnames = [os.path.join(caiman_datadir(), 'example_movies/Sue_2x_3000_40_-46.nwb')] # estimates save path can be same or different from raw data path - save_path = os.path.join(caiman_datadir(), 'example_movies/Sue_2x_3000_40_-46_CNMF_estimates.nwb') - # filename to be created or processed - # dataset dependent parameters - fr = 15. # imaging rate in frames per second - decay_time = 0.4 # length of a typical transient in seconds - - #%% load the file and save it in the NWB format (if it doesn't exist already) - if not os.path.exists(fnames[0]): - fnames_orig = 'Sue_2x_3000_40_-46.tif' # filename to be processed - if fnames_orig in ['Sue_2x_3000_40_-46.tif', 'demoMovie.tif']: - fnames_orig = [download_demo(fnames_orig)] - orig_movie = cm.load(fnames_orig, fr=fr) + save_path = os.path.splitext(fnames[0])[0] + '_CNMF_estimates.nwb' - # save file in NWB format with various additional info - orig_movie.save(fnames[0], sess_desc='test', identifier='demo 1', - imaging_plane_description='single plane', - emission_lambda=520.0, indicator='GCAMP6f', - location='parietal cortex', - experimenter='Sue Ann Koay', lab_name='Tank Lab', - institution='Princeton U', - experiment_description='Experiment Description', - session_id='Session 1', - var_name_hdf5='TwoPhotonSeries') - - #%% First setup some parameters for data and motion correction - # motion correction parameters - dxy = (2., 2.) # spatial resolution in x and y in (um per pixel) - # note the lower than usual spatial resolution here - max_shift_um = (12., 12.) # maximum shift in um - patch_motion_um = (100., 100.) # patch size for non-rigid correction in um - pw_rigid = True # flag to select rigid vs pw_rigid motion correction - # maximum allowed rigid shift in pixels - max_shifts = [int(a/b) for a, b in zip(max_shift_um, dxy)] - # start a new patch for pw-rigid motion correction every x pixels - strides = tuple([int(a/b) for a, b in zip(patch_motion_um, dxy)]) - # overlap between patches (size of patch in pixels: strides+overlaps) - overlaps = (24, 24) - # maximum deviation allowed for patch with respect to rigid shifts - max_deviation_rigid = 3 - - mc_dict = { - 'fnames': fnames, - 'fr': fr, - 'decay_time': decay_time, - 'dxy': dxy, - 'pw_rigid': pw_rigid, - 'max_shifts': max_shifts, - 'strides': strides, - 'overlaps': overlaps, - 'max_deviation_rigid': max_deviation_rigid, - 'border_nan': 'copy', - #'var_name_hdf5': 'acquisition/TwoPhotonSeries' - 'var_name_hdf5': 'TwoPhotonSeries' - } - opts = params.CNMFParams(params_dict=mc_dict) - - # %% play the movie (optional) + # We're done with all the file input parts + + m_orig = caiman.load_movie_chain(fnames, var_name_hdf5=opts.data['var_name_hdf5']) + + # play the movie (optional) # playing the movie using opencv. It requires loading the movie in memory. # To close the video press q - display_images = False - if display_images: - m_orig = cm.load_movie_chain(fnames, var_name_hdf5=opts.data['var_name_hdf5']) + + if not cfg.no_play: ds_ratio = 0.2 moviehandle = m_orig.resize(1, 1, ds_ratio) moviehandle.play(q_max=99.5, fr=60, magnification=2) - # %% start a cluster for parallel processing - c, dview, n_processes = cm.cluster.setup_cluster( - backend='multiprocessing', n_processes=None, single_thread=False) + # start a cluster for parallel processing + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) - # %%% MOTION CORRECTION + # Motion Correction # first we create a motion correction object with the specified parameters mc = MotionCorrect(fnames, dview=dview, var_name_hdf5=opts.data['var_name_hdf5'], **opts.get_group('motion')) # note that the file is not loaded in memory - # %% Run (piecewise-rigid motion) correction using NoRMCorre + # Run (piecewise-rigid motion) correction using NoRMCorre mc.motion_correct(save_movie=True) - # %% compare with original movie - if display_images: - m_orig = cm.load_movie_chain(fnames, var_name_hdf5=opts.data['var_name_hdf5']) - m_els = cm.load(mc.mmap_file) + # compare with original movie + if not cfg.no_play: + m_orig = caiman.load_movie_chain(fnames, var_name_hdf5=opts.data['var_name_hdf5']) + m_els = caiman.load(mc.mmap_file) ds_ratio = 0.2 - moviehandle = cm.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie, + moviehandle = caiman.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie, m_els.resize(1, 1, ds_ratio)], axis=2) moviehandle.play(fr=60, q_max=99.5, magnification=2) # press q to exit - # %% MEMORY MAPPING - border_to_0 = 0 if mc.border_nan == 'copy' else mc.border_to_0 + # MEMORY MAPPING + if mc.border_nan == 'copy': + border_to_0 = 0 + else: + border_to_0 = mc.border_to_0 # you can include the boundaries of the FOV if you used the 'copy' option # during motion correction, although be careful about the components near # the boundaries # memory map the file in order 'C' - fname_new = cm.save_memmap(mc.mmap_file, base_name='memmap_', order='C', + fname_new = caiman.save_memmap(mc.mmap_file, base_name='memmap_', order='C', border_to_0=border_to_0) # exclude borders # now load the file - Yr, dims, T = cm.load_memmap(fname_new) + Yr, dims, T = caiman.load_memmap(fname_new) images = np.reshape(Yr.T, [T] + list(dims), order='F') # load frames in python format (T x X x Y) - # %% restart cluster to clean up memory - cm.stop_server(dview=dview) - c, dview, n_processes = cm.cluster.setup_cluster( - backend='multiprocessing', n_processes=None, single_thread=False) - - # %% parameters for source extraction and deconvolution - p = 1 # order of the autoregressive system - gnb = 2 # number of global background components - merge_thr = 0.85 # merging threshold, max correlation allowed - rf = 15 - # half-size of the patches in pixels. e.g., if rf=25, patches are 50x50 - stride_cnmf = 6 # amount of overlap between the patches in pixels - K = 4 # number of components per patch - gSig = [4, 4] # expected half size of neurons in pixels - # initialization method (if analyzing dendritic data using 'sparse_nmf') - method_init = 'greedy_roi' - ssub = 2 # spatial subsampling during initialization - tsub = 2 # temporal subsampling during initialization - - # parameters for component evaluation - opts_dict = {'fnames': fnames, - 'fr': fr, - 'nb': gnb, - 'rf': rf, - 'K': K, - 'gSig': gSig, - 'stride': stride_cnmf, - 'method_init': method_init, - 'rolling_sum': True, - 'merge_thr': merge_thr, - 'n_processes': n_processes, - 'only_init': True, - 'ssub': ssub, - 'tsub': tsub} - - opts.change_params(params_dict=opts_dict); - - # %% RUN CNMF ON PATCHES + # restart cluster to clean up memory + caiman.stop_server(dview=dview) + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) + + # RUN CNMF ON PATCHES # First extract spatial and temporal components on patches and combine them # for this step deconvolution is turned off (p=0) - opts.change_params({'p': 0}) + saved_p = opts.preprocess['p'] # Save the deconvolution parameter for later restoration + opts.change_params({'preprocess': {'p': 0}, 'temporal': {'p': 0}}) cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) cnm = cnm.fit(images) - # %% ALTERNATE WAY TO RUN THE PIPELINE AT ONCE (optional) + # ALTERNATE WAY TO RUN THE PIPELINE AT ONCE (optional) # you can also perform the motion correction plus cnmf fitting steps # simultaneously after defining your parameters object using # cnm1 = cnmf.CNMF(n_processes, params=opts, dview=dview) # cnm1.fit_file(motion_correct=True) - # %% plot contours of found components - Cn = cm.local_correlations(images, swap_dim=False) + # plot contours of found components + Cn = caiman.summary_images.local_correlations(images, swap_dim=False) Cn[np.isnan(Cn)] = 0 - cnm.estimates.plot_contours(img=Cn) + if not cfg.no_play: + cnm.estimates.plot_contours(img=Cn) plt.title('Contour plots of found components') - #%% save results in a separate file (just for demonstration purposes) + # save results in a separate file (just for demonstration purposes) cnm.estimates.Cn = Cn - cnm.save(fname_new[:-4]+'hdf5') - #cm.movie(Cn).save(fname_new[:-5]+'_Cn.tif') + cnm.save(fname_new[:-4] + 'hdf5') # FIXME - # %% RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution - cnm.params.change_params({'p': p}) + # RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution + cnm.params.change_params({'preprocess': {'p': saved_p}, 'temporal': {'p': saved_p}}) # Restore deconvolution cnm2 = cnm.refit(images, dview=dview) - # %% COMPONENT EVALUATION - # the components are evaluated in three ways: - # a) the shape of each component must be correlated with the data - # b) a minimum peak SNR is required over the length of a transient - # c) each shape passes a CNN based classifier - min_SNR = 2 # signal to noise ratio for accepting a component - rval_thr = 0.85 # space correlation threshold for accepting a component - cnn_thr = 0.99 # threshold for CNN based classifier - cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected - - cnm2.params.set('quality', {'decay_time': decay_time, - 'min_SNR': min_SNR, - 'rval_thr': rval_thr, - 'use_cnn': True, - 'min_cnn_thr': cnn_thr, - 'cnn_lowest': cnn_lowest}); cnm2.estimates.evaluate_components(images, cnm2.params, dview=dview) - #%% Save new estimates object (after adding correlation image) + # Save new estimates object (after adding correlation image) cnm2.estimates.Cn = Cn cnm2.save(fname_new[:-4] + 'hdf5') - # %% PLOT COMPONENTS - cnm2.estimates.plot_contours(img=Cn, idx=cnm2.estimates.idx_components) + if not cfg.no_play: + # Plot Components + cnm2.estimates.plot_contours(img=Cn, idx=cnm2.estimates.idx_components) - # %% VIEW TRACES (accepted and rejected) - if display_images: + # View Traces (accepted and rejected) cnm2.estimates.view_components(images, img=Cn, idx=cnm2.estimates.idx_components) cnm2.estimates.view_components(images, img=Cn, idx=cnm2.estimates.idx_components_bad) - #%% update object with selected components (optional) + + # update object with selected components (optional) # cnm2.estimates.select_components(use_object=True) - #%% Extract DF/F values + # Extract DF/F values cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250) - #%% Show final traces - cnm2.estimates.view_components(img=Cn) + # Show final traces + if not cfg.no_play: + cnm2.estimates.view_components(img=Cn) - #%% reconstruct denoised movie (press q to exit) - if display_images: + # reconstruct denoised movie (press q to exit) + if not cfg.no_play: cnm2.estimates.play_movie(images, q_max=99.9, gain_res=2, magnification=2, bpx=border_to_0, include_bck=False) # background not shown - #%% STOP CLUSTER and clean up log files - cm.stop_server(dview=dview) - log_files = glob.glob('*_LOG_*') - for log_file in log_files: - os.remove(log_file) - - #%% save the results in the original NWB file + # Stop the cluster and clean up log files + caiman.stop_server(dview=dview) + + if not cfg.keep_logs: + log_files = glob.glob('*_LOG_*') + for log_file in log_files: + os.remove(log_file) + + # save the results in the original NWB file cnm2.estimates.save_NWB(save_path, imaging_rate=fr, session_start_time=datetime.now(tzlocal()), raw_data_file=fnames[0]) -# %% -# This is to mask the differences between running this demo in IDE -# versus from the CLI -if __name__ == "__main__": - main() +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate 2P Pipeline using batch algorithm, with nwb files") + parser.add_argument("--configfile", default=os.path.join(caiman.paths.caiman_datadir(), 'demos', 'general', 'params_demo_pipeline_NWB.json'), help="JSON Configfile for Caiman parameters") + parser.add_argument("--keep_logs", action="store_true", help="Keep temporary logfiles") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") + parser.add_argument("--cluster_nproc", type=int, default=None, help="Override automatic selection of number of workers to use") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() + +######## +main() + diff --git a/demos/general/demo_pipeline_cnmfE.py b/demos/general/demo_pipeline_cnmfE.py index e61178876..3b3d102bf 100755 --- a/demos/general/demo_pipeline_cnmfE.py +++ b/demos/general/demo_pipeline_cnmfE.py @@ -9,224 +9,116 @@ CNMF-E algorithm for source extraction (as opposed to plain CNMF). Check the companion paper for more details. -You can also run a large part of the pipeline with a single method -(cnmf.fit_file) See inside for details. - -Demo is also available as a jupyter notebook (see demo_pipeline_cnmfE.ipynb) +See also the jupyter notebook demo_pipeline_cnmfE.ipynb """ -#%% -from IPython import get_ipython + +import argparse +import glob import logging import matplotlib.pyplot as plt import numpy as np +import os -try: - if __IPYTHON__: - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: - pass - -import caiman as cm +import caiman +from caiman.motion_correction import MotionCorrect from caiman.source_extraction import cnmf +from caiman.source_extraction.cnmf import params as params from caiman.utils.utils import download_demo from caiman.utils.visualization import inspect_correlation_pnr -from caiman.motion_correction import MotionCorrect -from caiman.source_extraction.cnmf import params as params -#%% -# Set up the logger; change this if you like. -# You can log to a file using the filename parameter, or make the output more or less -# verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]"\ - "[%(process)d] %(message)s", - level=logging.WARNING) - # filename="/tmp/caiman.log" - -#%% -def main(): - pass # For compatibility between running under an IDE and the CLI - - #%% start the cluster - try: - cm.stop_server() # stop it if it was running - except(): - pass - c, dview, n_processes = cm.cluster.setup_cluster(backend='multiprocessing', - n_processes=24, # number of process to use, if you go out of memory try to reduce this one - single_thread=False) - #%% First setup some parameters for motion correction +def main(): + cfg = handle_args() + + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + opts = params.CNMFParams(params_from_file=cfg.configfile) + + if cfg.input is not None: + # If no input is specified, use sample data, downloading if necessary + opts.change_params({"data": {"fnames": cfg.input}}) + if not opts.data['fnames']: # Set neither by CLI arg nor through JSON, so use default data + fnames = [download_demo('data_endoscope.tif')] + opts.change_params({"data": {"fnames": fnames}}) + + # First set up some parameters for data and motion correction # dataset dependent parameters - fnames = ['data_endoscope.tif'] # filename to be processed - fnames = [download_demo(fnames[0])] # download file if not already present - filename_reorder = fnames - fr = 10 # movie frame rate - decay_time = 0.4 # length of a typical transient in seconds - - # motion correction parameters - motion_correct = True # flag for motion correction - pw_rigid = False # flag for pw-rigid motion correction - - gSig_filt = (3, 3) # size of filter, in general gSig (see below), - # change this one if algorithm does not work - max_shifts = (5, 5) # maximum allowed rigid shift - strides = (48, 48) # start a new patch for pw-rigid motion correction every x pixels - overlaps = (24, 24) # overlap between patches (size of patch strides+overlaps) - # maximum deviation allowed for patch with respect to rigid shifts - max_deviation_rigid = 3 - border_nan = 'copy' - - mc_dict = { - 'fnames': fnames, - 'fr': fr, - 'decay_time': decay_time, - 'pw_rigid': pw_rigid, - 'max_shifts': max_shifts, - 'gSig_filt': gSig_filt, - 'strides': strides, - 'overlaps': overlaps, - 'max_deviation_rigid': max_deviation_rigid, - 'border_nan': border_nan - } - - opts = params.CNMFParams(params_dict=mc_dict) - - # %% MOTION CORRECTION - # The pw_rigid flag set above, determines where to use rigid or pw-rigid + + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) + # Motion Correction + # The motion:pw_rigid parameter determines where to use rigid or pw-rigid # motion correction - if motion_correct: + if not cfg.no_motion_correct: # do motion correction rigid - mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion')) + mc = MotionCorrect(opts.data['fnames'], dview=dview, **opts.get_group('motion')) mc.motion_correct(save_movie=True) - fname_mc = mc.fname_tot_els if pw_rigid else mc.fname_tot_rig - if pw_rigid: + fname_mc = mc.fname_tot_els if opts.motion['pw_rigid'] else mc.fname_tot_rig + if opts.motion['pw_rigid']: bord_px = np.ceil(np.maximum(np.max(np.abs(mc.x_shifts_els)), np.max(np.abs(mc.y_shifts_els)))).astype(int) else: bord_px = np.ceil(np.max(np.abs(mc.shifts_rig))).astype(int) - plt.subplot(1, 2, 1); plt.imshow(mc.total_template_rig) # % plot template - plt.subplot(1, 2, 2); plt.plot(mc.shifts_rig) # % plot rigid shifts + plt.subplot(1, 2, 1); plt.imshow(mc.total_template_rig) # plot template + plt.subplot(1, 2, 2); plt.plot(mc.shifts_rig) # plot rigid shifts plt.legend(['x shifts', 'y shifts']) plt.xlabel('frames') plt.ylabel('pixels') - bord_px = 0 if border_nan == 'copy' else bord_px - fname_new = cm.save_memmap(fname_mc, base_name='memmap_', order='C', + bord_px = 0 if opts.motion['border_nan'] == 'copy' else bord_px + fname_new = caiman.save_memmap(fname_mc, base_name='memmap_', order='C', border_to_0=bord_px) else: # if no motion correction just memory map the file - fname_new = cm.save_memmap(filename_reorder, base_name='memmap_', + fname_new = caiman.save_memmap(opts.data['fnames'], base_name='memmap_', order='C', border_to_0=0, dview=dview) # load memory mappable file - Yr, dims, T = cm.load_memmap(fname_new) + Yr, dims, T = caiman.load_memmap(fname_new) images = Yr.T.reshape((T,) + dims, order='F') - # %% Parameters for source extraction and deconvolution (CNMF-E algorithm) - p = 1 # order of the autoregressive system - K = None # upper bound on number of components per patch, in general None for 1p data - gSig = (3, 3) # gaussian width of a 2D gaussian kernel, which approximates a neuron - gSiz = (13, 13) # average diameter of a neuron, in general 4*gSig+1 + # Parameters for source extraction and deconvolution (CNMF-E algorithm) Ain = None # possibility to seed with predetermined binary masks - merge_thr = .7 # merging threshold, max correlation allowed - rf = 40 # half-size of the patches in pixels. e.g., if rf=40, patches are 80x80 - stride_cnmf = 20 # amount of overlap between the patches in pixels - # (keep it at least large as gSiz, i.e 4 times the neuron size gSig) - tsub = 2 # downsampling factor in time for initialization, - # increase if you have memory problems - ssub = 1 # downsampling factor in space for initialization, - # increase if you have memory problems - # you can pass them here as boolean vectors - low_rank_background = None # None leaves background of each patch intact, - # True performs global low-rank approximation if gnb>0 - gnb = 0 # number of background components (rank) if positive, - # else exact ring model with following settings - # gnb= 0: Return background as b and W - # gnb=-1: Return full rank background B - # gnb<-1: Don't return background - nb_patch = 0 # number of background components (rank) per patch if gnb>0, - # else it is set automatically - min_corr = .8 # min peak value from correlation image - min_pnr = 10 # min peak to noise ration from PNR image - ssub_B = 2 # additional downsampling factor in space for background - ring_size_factor = 1.4 # radius of ring is gSiz*ring_size_factor - - opts.change_params(params_dict={'dims': dims, - 'method_init': 'corr_pnr', # use this for 1 photon - 'K': K, - 'gSig': gSig, - 'gSiz': gSiz, - 'merge_thr': merge_thr, - 'p': p, - 'tsub': tsub, - 'ssub': ssub, - 'rf': rf, - 'stride': stride_cnmf, - 'only_init': True, # set it to True to run CNMF-E - 'nb': gnb, - 'nb_patch': nb_patch, - 'method_deconvolution': 'oasis', # could use 'cvxpy' alternatively - 'low_rank_background': low_rank_background, - 'update_background_components': True, # sometimes setting to False improve the results - 'min_corr': min_corr, - 'min_pnr': min_pnr, - 'normalize_init': False, # just leave as is - 'center_psf': True, # leave as is for 1 photon - 'ssub_B': ssub_B, - 'ring_size_factor': ring_size_factor, - 'del_duplicates': True, # whether to remove duplicates from initialization + + opts.change_params(params_dict={'dims': dims, # we rework the source files 'border_pix': bord_px}) # number of pixels to not consider in the borders) - # %% compute some summary images (correlation and peak to noise) + # compute some summary images (correlation and peak to noise) # change swap dim if output looks weird, it is a problem with tiffile - cn_filter, pnr = cm.summary_images.correlation_pnr(images[::1], gSig=gSig[0], swap_dim=False) + cn_filter, pnr = caiman.summary_images.correlation_pnr(images[::1], gSig=opts.init['gSig'][0], swap_dim=False) # if your images file is too long this computation will take unnecessarily # long time and consume a lot of memory. Consider changing images[::1] to # images[::5] or something similar to compute on a subset of the data # inspect the summary images and set the parameters inspect_correlation_pnr(cn_filter, pnr) - # print parameters set above, modify them if necessary based on summary images - print(min_corr) # min correlation of peak (from correlation image) - print(min_pnr) # min peak to noise ratio + print(f"Minimum correlation: {opts.init['min_corr']}") # min correlation of peak (from correlation image) + print(f"Minimum peak to noise ratio: {opts.init['min_pnr']}") # min peak to noise ratio - # %% RUN CNMF ON PATCHES + # Run CMNF in patches cnm = cnmf.CNMF(n_processes=n_processes, dview=dview, Ain=Ain, params=opts) cnm.fit(images) - # %% ALTERNATE WAY TO RUN THE PIPELINE AT ONCE (optional -- commented out) - # you can also perform the motion correction plus cnmf fitting steps - # simultaneously after defining your parameters object using - # cnm1 = cnmf.CNMF(n_processes, params=opts, dview=dview) - # cnm1.fit_file(motion_correct=True) - - # %% Quality Control: DISCARD LOW QUALITY COMPONENTS - min_SNR = 2.5 # adaptive way to set threshold on the transient size - r_values_min = 0.85 # threshold on space consistency (if you lower more components - # will be accepted, potentially with worst quality) - cnm.params.set('quality', {'min_SNR': min_SNR, - 'rval_thr': r_values_min, - 'use_cnn': False}) + # Quality Control: DISCARD LOW QUALITY COMPONENTS cnm.estimates.evaluate_components(images, cnm.params, dview=dview) print(' ***** ') - print('Number of total components: ', len(cnm.estimates.C)) - print('Number of accepted components: ', len(cnm.estimates.idx_components)) + print(f"Number of total components: {len(cnm.estimates.C)}") + print(f"Number of accepted components: {len(cnm.estimates.idx_components)}") - # %% PLOT COMPONENTS - cnm.dims = dims - display_images = True # Set to true to show movies and images - if display_images: + # Play result movies + if not cfg.no_play: + cnm.dims = dims cnm.estimates.plot_contours(img=cn_filter, idx=cnm.estimates.idx_components) cnm.estimates.view_components(images, idx=cnm.estimates.idx_components) - - # %% MOVIES - display_images = False # Set to true to show movies and images - if display_images: # fully reconstructed movie cnm.estimates.play_movie(images, q_max=99.5, magnification=2, include_bck=True, gain_res=10, bpx=bord_px) @@ -234,11 +126,26 @@ def main(): cnm.estimates.play_movie(images, q_max=99.9, magnification=2, include_bck=False, gain_res=4, bpx=bord_px) - # %% STOP SERVER - cm.stop_server(dview=dview) - -# %% This is to mask the differences between running this demo in IDE -# versus from the CLI -if __name__ == "__main__": - main() + # Stop the cluster and clean up log files + caiman.stop_server(dview=dview) + + if not cfg.keep_logs: + log_files = glob.glob('*_LOG_*') + for log_file in log_files: + os.remove(log_file) + +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate CNMFE Pipeline") + parser.add_argument("--configfile", default=os.path.join(caiman.paths.caiman_datadir(), 'demos', 'general', 'params_demo_pipeline_cnmfE.json'), help="JSON Configfile for Caiman parameters") + parser.add_argument("--no_motion_correct", action='store_true', help="Set to disable motion correction") + parser.add_argument("--keep_logs", action="store_true", help="Keep temporary logfiles") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") + parser.add_argument("--cluster_nproc", type=int, default=None, help="Override automatic selection of number of workers to use") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--logfile", help="If specified, log to the named file") + return parser.parse_args() + +######## +main() diff --git a/demos/general/demo_pipeline_voltage_imaging.py b/demos/general/demo_pipeline_voltage_imaging.py old mode 100644 new mode 100755 index daa9f254e..20ad228e1 --- a/demos/general/demo_pipeline_voltage_imaging.py +++ b/demos/general/demo_pipeline_voltage_imaging.py @@ -1,4 +1,5 @@ #!/usr/bin/env python + """ Demo pipeline for processing voltage imaging data. The processing pipeline includes motion correction, memory mapping, segmentation, denoising and source @@ -7,11 +8,11 @@ Dataset courtesy of Karel Svoboda Lab (Janelia Research Campus). author: @caichangjia """ -#%% + +import argparse import cv2 import glob import h5py -from IPython import get_ipython import logging import matplotlib.pyplot as plt import numpy as np @@ -22,45 +23,41 @@ except: pass -try: - if __IPYTHON__: - ipython = get_ipython() - ipython.run_line_magic('load_ext', 'autoreload') - ipython.run_line_magic('autoreload', '2') - ipython.run_line_magic('matplotlib', 'qt') -except NameError: - pass - -import caiman as cm +import caiman from caiman.motion_correction import MotionCorrect from caiman.paths import caiman_datadir from caiman.source_extraction.volpy import utils from caiman.source_extraction.volpy.volparams import volparams from caiman.source_extraction.volpy.volpy import VOLPY -from caiman.summary_images import local_correlations_movie_offline -from caiman.summary_images import mean_image +from caiman.summary_images import local_correlations_movie_offline, mean_image from caiman.utils.utils import download_demo, download_model -# %% -# Set up the logger (optional); change this if you like. -# You can log to a file using the filename parameter, or make the output more -# or less verbose by setting level to logging.DEBUG, logging.INFO, -# logging.WARNING, or logging.ERROR -logging.basicConfig(format= - "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s]" \ - "[%(process)d] %(message)s", - level=logging.INFO) -# %% def main(): - pass # For compatibility between running using an IDE and the CLI + cfg = handle_args() - # %% Load demo movie and ROIs - fnames = download_demo('demo_voltage_imaging.hdf5', 'volpy') # file path to movie file (will download if not present) - path_ROIs = download_demo('demo_voltage_imaging_ROIs.hdf5', 'volpy') # file path to ROIs file (will download if not present) + if cfg.logfile: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO, + filename=cfg.logfile) + # You can make the output more or less verbose by setting level to logging.DEBUG, logging.INFO, logging.WARNING, or logging.ERROR + else: + logging.basicConfig(format= + "[%(filename)s:%(funcName)20s():%(lineno)s] %(message)s", + level=logging.INFO) + + if cfg.input is None: + # If no input is specified, use sample data, downloading if necessary + fnames = [download_demo('demo_voltage_imaging.hdf5', 'volpy')] # XXX do we need to use a separate directory? + path_ROIs = [download_demo('demo_voltage_imaging_ROIs.hdf5', 'volpy')] + else: + fnames = cfg.input + path_ROIs = cfg.pathinput + + # Load demo movie and ROIs file_dir = os.path.split(fnames)[0] - #%% dataset dependent parameters # dataset dependent parameters fr = 400 # sample rate of the movie @@ -74,33 +71,29 @@ def main(): max_deviation_rigid = 3 # maximum deviation allowed for patch with respect to rigid shifts border_nan = 'copy' - opts_dict = { - 'fnames': fnames, - 'fr': fr, - 'pw_rigid': pw_rigid, - 'max_shifts': max_shifts, - 'gSig_filt': gSig_filt, - 'strides': strides, - 'overlaps': overlaps, - 'max_deviation_rigid': max_deviation_rigid, - 'border_nan': border_nan - } + params_dict = {'fnames': fnames, + 'fr': fr, + 'pw_rigid': pw_rigid, + 'max_shifts': max_shifts, + 'gSig_filt': gSig_filt, + 'strides': strides, + 'overlaps': overlaps, + 'max_deviation_rigid': max_deviation_rigid, + 'border_nan': border_nan} - opts = volparams(params_dict=opts_dict) + opts = volparams(params_dict=params_dict) - # %% play the movie (optional) + # play the movie (optional) # playing the movie using opencv. It requires loading the movie in memory. # To close the movie press q - display_images = False - if display_images: - m_orig = cm.load(fnames) + if not cfg.no_play: + m_orig = caiman.load(fnames) ds_ratio = 0.2 moviehandle = m_orig.resize(1, 1, ds_ratio) moviehandle.play(q_max=99.5, fr=40, magnification=4) - # %% start a cluster for parallel processing - c, dview, n_processes = cm.cluster.setup_cluster( - backend='multiprocessing', n_processes=None, single_thread=False) + # start a cluster for parallel processing + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) # %%% MOTION CORRECTION # first we create a motion correction object with the specified parameters @@ -115,16 +108,16 @@ def main(): mc.mmap_file = [os.path.join(file_dir, mc_list[0])] print(f'reuse previously saved motion corrected file:{mc.mmap_file}') - # %% compare with original movie - if display_images: - m_orig = cm.load(fnames) - m_rig = cm.load(mc.mmap_file) + # compare with original movie + if not cfg.no_play: + m_orig = caiman.load(fnames) + m_rig = caiman.load(mc.mmap_file) ds_ratio = 0.2 - moviehandle = cm.concatenate([m_orig.resize(1, 1, ds_ratio), + moviehandle = caiman.concatenate([m_orig.resize(1, 1, ds_ratio), m_rig.resize(1, 1, ds_ratio)], axis=2) moviehandle.play(fr=40, q_max=99.5, magnification=4) # press q to exit - # %% MEMORY MAPPING + # MEMORY MAPPING do_memory_mapping = True if do_memory_mapping: border_to_0 = 0 if mc.border_nan == 'copy' else mc.border_to_0 @@ -133,15 +126,15 @@ def main(): # the boundaries # memory map the file in order 'C' - fname_new = cm.save_memmap_join(mc.mmap_file, base_name='memmap_' + os.path.splitext(os.path.split(fnames)[-1])[0], + fname_new = caiman.save_memmap_join(mc.mmap_file, base_name='memmap_' + os.path.splitext(os.path.split(fnames)[-1])[0], add_to_mov=border_to_0, dview=dview) # exclude border else: mmap_list = [file for file in os.listdir(file_dir) if ('memmap_' + os.path.splitext(os.path.split(fnames)[-1])[0]) in file] fname_new = os.path.join(file_dir, mmap_list[0]) - print(f'reuse previously saved memory mapping file:{fname_new}') + print(f'reuse previously saved memory mapping file: {fname_new}') - # %% SEGMENTATION + # SEGMENTATION # create summary images img = mean_image(mc.mmap_file[0], window = 1000, dview=dview) img = (img-np.mean(img))/np.std(img) @@ -154,30 +147,23 @@ def main(): img_corr = (Cn-np.mean(Cn))/np.std(Cn) summary_images = np.stack([img, img, img_corr], axis=0).astype(np.float32) # save summary images which are used in the VolPy GUI - cm.movie(summary_images).save(fnames[:-5] + '_summary_images.tif') + caiman.movie(summary_images).save(fnames[:-5] + '_summary_images.tif') fig, axs = plt.subplots(1, 2) axs[0].imshow(summary_images[0]); axs[1].imshow(summary_images[2]) axs[0].set_title('mean image'); axs[1].set_title('corr image'); - #%% methods for segmentation - methods_list = ['manual_annotation', # manual annotations need prepared annotated datasets in the same format as demo_voltage_imaging_ROIs.hdf5 - 'maskrcnn', # Mask R-CNN is a convolutional neural network trained for detecting neurons in summary images - 'gui_annotation'] # use VolPy GUI to correct outputs of Mask R-CNN or annotate new datasets - - method = methods_list[0] - if method == 'manual_annotation': + if cfg.method == 'manual_annotation': with h5py.File(path_ROIs, 'r') as fl: ROIs = fl['mov'][()] - elif method == 'maskrcnn': + elif cfg.method == 'maskrcnn': weights_path = download_model('mask_rcnn') ROIs = utils.mrcnn_inference(img=summary_images.transpose([1, 2, 0]), size_range=[5, 22], weights_path=weights_path, display_result=True) # size parameter decides size range of masks to be selected - cm.movie(ROIs).save(fnames[:-5] + '_mrcnn_ROIs.hdf5') + caiman.movie(ROIs).save(fnames[:-5] + '_mrcnn_ROIs.hdf5') - elif method == 'gui_annotation': + elif cfg.method == 'gui_annotation': # run volpy_gui.py file in the caiman/source_extraction/volpy folder - # or run the following in the ipython: %run volpy_gui.py gui_ROIs = caiman_datadir() + '/example_movies/volpy/gui_roi.hdf5' with h5py.File(gui_ROIs, 'r') as fl: ROIs = fl['mov'][()] @@ -186,12 +172,12 @@ def main(): axs[0].imshow(summary_images[0]); axs[1].imshow(ROIs.sum(0)) axs[0].set_title('mean image'); axs[1].set_title('masks') - # %% restart cluster to clean up memory - cm.stop_server(dview=dview) - c, dview, n_processes = cm.cluster.setup_cluster( - backend='multiprocessing', n_processes=None, single_thread=False, maxtasksperchild=1) + # restart cluster to clean up memory + caiman.stop_server(dview=dview) + c, dview, n_processes = caiman.cluster.setup_cluster(backend=cfg.cluster_backend, n_processes=cfg.cluster_nproc) + - # %% parameters for trace denoising and spike extraction + # parameters for trace denoising and spike extraction ROIs = ROIs # region of interests index = list(range(len(ROIs))) # index of neurons weights = None # if None, use ROIs for initialization; to reuse weights check reuse weights block @@ -232,27 +218,26 @@ def main(): 'weight_update': weight_update, 'n_iter': n_iter} - opts.change_params(params_dict=opts_dict); + opts.change_params(params_dict=params_dict); - #%% TRACE DENOISING AND SPIKE DETECTION + # TRACE DENOISING AND SPIKE DETECTION vpy = VOLPY(n_processes=n_processes, dview=dview, params=opts) vpy.fit(n_processes=n_processes, dview=dview); - #%% visualization - display_images = True - if display_images: + # visualization + if not cfg.no_play: print(np.where(vpy.estimates['locality'])[0]) # neurons that pass locality test idx = np.where(vpy.estimates['locality'] > 0)[0] utils.view_components(vpy.estimates, img_corr, idx) - #%% reconstructed movie + # reconstructed movie # note the negative spatial weights are cutoff - if display_images: + if not cfg.no_play: mv_all = utils.reconstructed_movie(vpy.estimates.copy(), fnames=mc.mmap_file, idx=idx, scope=(0,1000), flip_signal=flip_signal) mv_all.play(fr=40, magnification=3) - #%% save the result in .npy format + # save the result in .npy format save_result = True if save_result: vpy.estimates['ROIs'] = ROIs @@ -260,7 +245,7 @@ def main(): save_name = f'volpy_{os.path.split(fnames)[1][:-5]}_{threshold_method}' np.save(os.path.join(file_dir, save_name), vpy.estimates) - #%% reuse weights + # reuse weights # set weights = reuse_weights in opts_dict dictionary estimates = np.load(os.path.join(file_dir, save_name+'.npy'), allow_pickle=True).item() reuse_weights = [] @@ -270,15 +255,26 @@ def main(): plt.figure(); plt.imshow(w);plt.colorbar(); plt.show() reuse_weights.append(w) - # %% STOP CLUSTER and clean up log files - cm.stop_server(dview=dview) - log_files = glob.glob('*_LOG_*') - for log_file in log_files: - os.remove(log_file) + # Stop the cluster and clean up log files + caiman.stop_server(dview=dview) + + if not cfg.keep_logs: + log_files = glob.glob('*_LOG_*') + for log_file in log_files: + os.remove(log_file) + +def handle_args(): + parser = argparse.ArgumentParser(description="Demonstrate voltage imaging pipeline") + parser.add_argument("--keep_logs", action="store_true", help="Keep temporary logfiles") + parser.add_argument("--no_play", action="store_true", help="Do not display results") + parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") + parser.add_argument("--cluster_nproc", type=int, default=None, help="Override automatic selection of number of workers to use") + parser.add_argument("--input", action="append", help="File(s) to work on, provide multiple times for more files") + parser.add_argument("--pathinput", action="append", help="Path input file(s) to work on, provide multiple times for more files") + parser.add_argument("--logfile", help="If specified, log to the named file") + parser.add_argument("--method", default="manual_annotation", choices=["manual_annotation", "maskrcnn", "gui_annotation"], help="Segmentation method") + return parser.parse_args() -# %% -# This is to mask the differences between running this demo in an IDE -# versus from the CLI -if __name__ == "__main__": - main() +######## +main() diff --git a/demos/general/params_demo_OnACID.json b/demos/general/params_demo_OnACID.json new file mode 100644 index 000000000..d871c6f98 --- /dev/null +++ b/demos/general/params_demo_OnACID.json @@ -0,0 +1,25 @@ +{ +"data": { + "decay_time": 0.75, + "fr": 10 + }, +"init": { + "gSig": [6, 6], + "K": 4, + "nb": 2 + }, +"online": { + "init_batch": 400, + "init_method": "cnmf", + "min_SNR": 1, + "sniper_mode": true, + "thresh_CNN_noisy": 0.65 + }, +"patch":{ + "rf": 16, + "stride": 3 + }, +"preprocess": { + "p": 1 + } +} diff --git a/demos/general/params_demo_OnACID_mesoscope.json b/demos/general/params_demo_OnACID_mesoscope.json new file mode 100644 index 000000000..623501a18 --- /dev/null +++ b/demos/general/params_demo_OnACID_mesoscope.json @@ -0,0 +1,41 @@ +{ + "data": { + "fr": 15, + "decay_time": 0.5 + }, + "init": { + "gSig": [3, 3], + "nb": 2, + "K": 2 + }, + "preprocess": { + "p": 1 + }, + "temporal": { + "p": 1 + }, + "online": { + "min_SNR": 1, + "rval_thr": 0.9, + "ds_factor": 1, + "motion_correct": true, + "init_batch": 200, + "init_method": "bare", + "normalize": true, + "sniper_mode": true, + "epochs": 1, + "max_shifts_online": 10, + "dist_shape_update": true, + "min_num_trial": 10 + }, + "motion": { + "pw_rigid": false + }, + "quality": { + "min_SNR": 2, + "rval_thr": 0.85, + "use_cnn": true, + "min_cnn_thr": 0.99, + "cnn_lowest": 0.1 + } +} diff --git a/demos/general/params_demo_caiman_basic_with_patches.json b/demos/general/params_demo_caiman_basic_with_patches.json new file mode 100644 index 000000000..7cd7bae8e --- /dev/null +++ b/demos/general/params_demo_caiman_basic_with_patches.json @@ -0,0 +1,31 @@ +{ +"data": { + "fr": 10, + "decay_time": 5.0 + }, +"init": { + "gSig": [6, 6], + "K": 4, + "nb": 2 + }, +"patch": { + "rf": 10, + "stride": 4 + }, +"merging": { + "merge_thr": 0.80 + }, +"preprocess": { + "p": 2 + }, +"temporal": { + "p": 2 + }, +"quality": { + "min_SNR": 2, + "rval_thr": 0.85, + "use_cnn": true, + "min_cnn_thr": 0.99, + "cnn_lowest": 0.1 + } +} diff --git a/demos/general/params_demo_caiman_basic_without_patches.json b/demos/general/params_demo_caiman_basic_without_patches.json new file mode 100644 index 000000000..229fffe0e --- /dev/null +++ b/demos/general/params_demo_caiman_basic_without_patches.json @@ -0,0 +1,27 @@ +{ +"data": { + "fr": 10, + "decay_time": 5.0 + }, +"init": { + "gSig": [6, 6], + "K": 30, + "nb": 2 + }, +"merging": { + "merge_thr": 0.80 + }, +"preprocess": { + "p": 2 + }, +"temporal": { + "p": 2 + }, +"quality": { + "min_SNR": 2, + "rval_thr": 0.85, + "use_cnn": true, + "min_cnn_thr": 0.99, + "cnn_lowest": 0.1 + } +} diff --git a/demos/general/params_demo_pipeline.json b/demos/general/params_demo_pipeline.json new file mode 100644 index 000000000..567f1d255 --- /dev/null +++ b/demos/general/params_demo_pipeline.json @@ -0,0 +1,45 @@ +{ + "data": { + "fr": 30, + "decay_time": 0.4, + "dxy": [2.0, 2.0] + }, + "init": { + "nb": 2, + "K": 4, + "gSig": [4, 4], + "method_init": "greedy_roi", + "rolling_sum": true, + "ssub": 2, + "tsub": 2 + }, + "motion": { + "pw_rigid": true, + "max_shifts": [6, 6], + "strides": [50, 50], + "overlaps": [24, 24], + "max_deviation_rigid": 3, + "border_nan": "copy" + }, + "preprocess": { + "p": 1 + }, + "temporal": { + "p": 1 + }, + "patch": { + "rf": 15, + "stride": 6, + "only_init": true + }, + "merging": { + "merge_thr": 0.85 + }, + "quality": { + "min_SNR": 2, + "rval_thr": 0.85, + "use_cnn": true, + "min_cnn_thr": 0.99, + "cnn_lowest": 0.1 + } +} diff --git a/demos/general/params_demo_pipeline_NWB.json b/demos/general/params_demo_pipeline_NWB.json new file mode 100644 index 000000000..6cca1d199 --- /dev/null +++ b/demos/general/params_demo_pipeline_NWB.json @@ -0,0 +1,45 @@ +{ + "data": { + "fr": 15.0, + "decay_time": 0.4, + "dxy": [2.0, 2.0] + }, + "init": { + "nb": 2, + "K": 4, + "gSig": [4, 4], + "method_init": "greedy_roi", + "rolling_sum": true, + "ssub": 2, + "tsub": 2 + }, + "motion": { + "pw_rigid": true, + "max_shifts": [6, 6], + "strides": [50, 50], + "overlaps": [24, 24], + "max_deviation_rigid": 3, + "border_nan": "copy" + }, + "preprocess": { + "p": 1 + }, + "temporal": { + "p": 1 + }, + "patch": { + "rf": 15, + "stride": 6, + "only_init": true + }, + "merging": { + "merge_thr": 0.85 + }, + "quality": { + "min_SNR": 2, + "rval_thr": 0.85, + "use_cnn": true, + "min_cnn_thr": 0.99, + "cnn_lowest": 0.1 + } +} diff --git a/demos/general/params_demo_pipeline_cnmfE.json b/demos/general/params_demo_pipeline_cnmfE.json new file mode 100644 index 000000000..8f68a2cea --- /dev/null +++ b/demos/general/params_demo_pipeline_cnmfE.json @@ -0,0 +1,56 @@ +{ + "data": { + "fr": 10, + "decay_time": 0.4 + }, + "init": { + "method_init": "corr_pnr", + "K": null, + "gSig": [3, 3], + "gSiz": [13, 13], + "ssub": 1, + "tsub": 2, + "nb": 0, + "min_corr": 0.8, + "min_pnr": 10, + "normalize_init": false, + "center_psf": true, + "ssub_B": 2, + "ring_size_factor": 1.4 + }, + "motion": { + "pw_rigid": false, + "max_shifts": [5, 5], + "gSig_filt": [3, 3], + "strides": [48, 48], + "overlaps": [24, 24], + "max_deviation_rigid": 3, + "border_nan": "copy" + }, + "preprocess": { + "p": 1 + }, + "temporal": { + "p": 1, + "method_deconvolution": "oasis" + }, + "spatial": { + "update_background_components": true + }, + "patch": { + "rf": 40, + "stride": 20, + "only_init": true, + "nb_patch": 0, + "low_rank_background": null, + "del_duplicates": true + }, + "merging": { + "merge_thr": 0.7 + }, + "quality": { + "min_SNR": 2.5, + "rval_thr": 0.85, + "use_cnn": false + } +} diff --git a/docs/source/file_formats.md b/docs/source/file_formats.md index 6e1a94246..4eeef523b 100644 --- a/docs/source/file_formats.md +++ b/docs/source/file_formats.md @@ -15,7 +15,7 @@ Notes for specific formats follow: Matlab (\*.mat) --------------- -Caiman may be able to handle Matlab's mat files, although it is not a desirable format. We rely on the scipy.io.matlab functions to handle files of this type, except for certain versions of the format which are actually hdf5 files +Caiman can handle Matlab's V2 mat files (which are HDF5 files in a certain format). If you need old-style support, reach out and we can consider adding it. TIFF (\*.tiff, \*.tif, \*.btf) ------------------------------ diff --git a/environment-minimal.yml b/environment-minimal.yml index 5f16cbd69..6ff26de47 100644 --- a/environment-minimal.yml +++ b/environment-minimal.yml @@ -1,19 +1,23 @@ channels: - conda-forge +- defaults dependencies: - python <=3.12 - av - cython -- future -- h5py +- h5py >=3.4.0 - holoviews >=1.16.2 +- ipython - ipyparallel +- ipywidgets - matplotlib - moviepy -- numpy <2.0.0 +- nose +- numpy <2.0.0,>=1.26 - numpydoc - opencv - peakutils +- pims - psutil - pynwb - scikit-image >=0.19.0 diff --git a/environment.yml b/environment.yml index 4f087c33d..6ae3e04e5 100644 --- a/environment.yml +++ b/environment.yml @@ -1,12 +1,13 @@ channels: - conda-forge +- defaults dependencies: - python <=3.12 - av - bokeh >=3.1.1 - coverage - cython -- h5py +- h5py >=3.4.0 - holoviews >=1.16.2 - ipykernel - ipython @@ -17,7 +18,7 @@ dependencies: - moviepy - mypy - nose -- numpy <2.0.0 +- numpy <2.0.0,>=1.26 - numpydoc - opencv - panel >=1.0.2