diff --git a/Dockerfile b/Dockerfile index 473576ad4..259f9b5e2 100644 --- a/Dockerfile +++ b/Dockerfile @@ -2,7 +2,7 @@ FROM continuumio/anaconda3 RUN conda config --set always_yes yes RUN conda update --yes conda -RUN apt-get install -y gcc g++ libgl1 +RUN apt-get update && apt-get install -y gcc g++ libgl1 RUN mkdir src && cd src && git clone -b dev https://github.com/flatironinstitute/CaImAn.git && cd CaImAn && conda env create -n caiman -f environment.yml && conda install --override-channels -c conda-forge -n caiman pip RUN /bin/bash -c "cd src/CaImAn && source activate caiman && /opt/conda/envs/caiman/bin/pip install ." RUN /bin/bash -c "source activate caiman && caimanmanager.py install" diff --git a/INSTALL-windows.md b/INSTALL-windows.md index 9936b6b79..95745d9d3 100755 --- a/INSTALL-windows.md +++ b/INSTALL-windows.md @@ -3,6 +3,7 @@ The Windows installation process differs more widely from installation on Linux ### Process * Increase the maximum size of your pagefile to 64G or more (http://www.tomshardware.com/faq/id-2864547/manage-virtual-memory-pagefile-windows.html ) - The Windows memmap interface is sensitive to the maximum setting and leaving it at the default can cause errors when processing larger datasets + * Remove any associations you may have made between .py files and an existing python interpreter or editor * Download and install Anaconda (Python 3.x, not 2.x) . Allow the installer to modify your PATH variable * Use Conda to install git (With "conda install git") - use of another commandline git is acceptable, but may lead to issues depending on default settings * Install Microsoft Build Tools for Visual Studio 2017 . Check the "Build Tools" box, and in the detailed view on the right check the "C/C++ CLI Tools" component too. The specifics of this occasionally change as Microsoft changes its products and website; you may need to go off-script. diff --git a/caiman/base/movies.py b/caiman/base/movies.py index 40489d170..8eac3293f 100644 --- a/caiman/base/movies.py +++ b/caiman/base/movies.py @@ -28,6 +28,7 @@ from matplotlib import animation import numpy as np import os +from PIL import Image # $ pip install pillow import pylab as pl import scipy.ndimage import scipy @@ -41,11 +42,11 @@ import sys import tifffile from tqdm import tqdm +from typing import List, Tuple import warnings from zipfile import ZipFile -from PIL import Image # $ pip install pillow -import caiman as cm +import caiman as cm from . import timeseries @@ -106,12 +107,6 @@ def __new__(cls, input_arr, **kwargs): else: raise Exception('Input must be an ndarray, use load instead!') - def motion_correction_online(self, max_shift_w=25, max_shift_h=25, init_frames_template=100, - show_movie=False, bilateral_blur=False, template=None, min_count=1000): - return motion_correct_online(self, max_shift_w=max_shift_w, max_shift_h=max_shift_h, - init_frames_template=init_frames_template, show_movie=show_movie, - bilateral_blur=bilateral_blur, template=template, min_count=min_count) - def apply_shifts_online(self, xy_shifts, save_base_name=None): # todo: todocument @@ -427,7 +422,7 @@ def linf(x, a, b): return a * x + b try: - p0 = (y[0] - y[-1], 1e-6, y[-1]) + p0:Tuple = (y[0] - y[-1], 1e-6, y[-1]) popt, _ = scipy.optimize.curve_fit(expf, x, y, p0=p0) y_fit = expf(x, *popt) except: @@ -825,7 +820,7 @@ def resize(self, fx=1, fy=1, fz=1, interpolation=cv2.INTER_AREA): max_els = 2**31 - 1 if elm > max_els: chunk_size = old_div((max_els), d) - new_m = [] + new_m:List = [] logging.debug('Resizing in chunks because of opencv bug') for chunk in range(0, T, chunk_size): logging.debug([chunk, np.minimum(chunk + chunk_size, T)]) @@ -1094,34 +1089,34 @@ def animate(i): for i in range(10): cv2.waitKey(100) -def load(file_name,fr=30,start_time=0,meta_data=None,subindices=None,shape=None, - var_name_hdf5 = 'mov', in_memory = False, is_behavior = False, bottom=0, - top=0, left=0, right=0, channel = None, outtype=np.float32): +def load(file_name, fr=30, start_time=0, meta_data=None, subindices=None, + shape=None, var_name_hdf5='mov', in_memory=False, is_behavior=False, + bottom=0, top=0, left=0, right=0, channel = None, outtype=np.float32): """ load movie from file. SUpports a variety of formats. tif, hdf5, npy and memory mapped. Matlab is experimental. Args: file_name: string name of file. Possible extensions are tif, avi, npy, (npz and hdf5 are usable only if saved by calblitz) - + fr: float frame rate - + start_time: float initial time for frame 1 - + meta_data: dict dictionary containing meta information about the movie - + subindices: iterable indexes for loading only portion of the movie - + shape: tuple of two values dimension of the movie along x and y if loading from a two dimensional numpy array - + num_frames_sub_idx: when reading sbx format (experimental and unstable) - + var_name_hdf5: str if loading from hdf5 name of the variable to load @@ -1140,39 +1135,45 @@ def load(file_name,fr=30,start_time=0,meta_data=None,subindices=None,shape=None, Exception 'File not found!' """ # case we load movie from file + if max(top, bottom, left, right) > 0 and type(file_name) is str: + file_name = [file_name] + if type(file_name) is list: if shape is not None: - raise Exception('shape not supported for multiple movie input') + logging.error('shape not supported for multiple movie input') return load_movie_chain(file_name,fr=fr, start_time=start_time, meta_data=meta_data, subindices=subindices, bottom=bottom, top=top, left=left, right=right, channel = channel, outtype=outtype) - if bottom != 0: - raise Exception('top bottom etc... not supported for single movie input') + if max(top, bottom, left, right) > 0: + logging.error('top bottom etc... not supported for single movie input') if channel is not None: - raise Exception('channel not supported for single movie input') + logging.error('channel not supported for single movie input') if os.path.exists(file_name): _, extension = os.path.splitext(file_name)[:2] extension = extension.lower() if extension == '.tif' or extension == '.tiff': # load avi file with tifffile.TiffFile(file_name) as tffl: + multi_page = True if tffl.series[0].shape[0] > 1 else False + if len(tffl.pages) == 1: + logging.warning('Your tif file is saved a single page' + + 'file. Performance will be affected') + multi_page = False if subindices is not None: if type(subindices) is list: - try: + if multi_page: input_arr = tffl.asarray(key=subindices[0])[:, subindices[1], subindices[2]] - except: - logging.warning('Your tif file is saved a single page file. Performance will be affected') + else: input_arr = tffl.asarray() input_arr = input_arr[subindices[0], subindices[1], subindices[2]] else: - try: + if multi_page: input_arr = tffl.asarray(key=subindices) - except: - logging.warning('Your tif file is saved a single page file. Performance will be affected') + else: input_arr = tffl.asarray() input_arr = input_arr[subindices] @@ -1287,25 +1288,7 @@ def rgb2gray(rgb): with np.load(file_name) as f: return movie(**f).astype(outtype) -# elif extension in ('.hdf5', '.h5'): -# with h5py.File(file_name, "r") as f: -# attrs = dict(f[var_name_hdf5].attrs) -# if meta_data in attrs: -# attrs['meta_data'] = cpk.loads(attrs['meta_data']) -# -# if subindices is None: -# return movie(f[var_name_hdf5], **attrs).astype(outtype) -# else: -# return movie(f[var_name_hdf5][subindices], **attrs).astype(outtype) - - elif extension == '.h5_at': - with h5py.File(file_name, "r") as f: - if subindices is None: - return movie(f['quietBlock'], fr=fr).astype(outtype) - else: - return movie(f['quietBlock'][subindices], fr=fr).astype(outtype) - - elif extension in ('.hdf5', '.h5'): + elif extension in ('.hdf5', '.h5', '.nwb'): if is_behavior: with h5py.File(file_name, "r") as f: kk = list(f.keys()) @@ -1322,22 +1305,22 @@ def rgb2gray(rgb): fkeys = list(f.keys()) if len(fkeys) == 1: var_name_hdf5 = fkeys[0] - if var_name_hdf5 in fkeys: + if var_name_hdf5 in f: if subindices is None: images = np.array(f[var_name_hdf5]).squeeze() - if images.ndim > 3: - images = images[:, 0] + #if images.ndim > 3: + # images = images[:, 0] else: images = np.array( f[var_name_hdf5][subindices]).squeeze() - if images.ndim > 3: - images = images[:, 0] + #if images.ndim > 3: + # images = images[:, 0] #input_arr = images return movie(images.astype(outtype)) else: logging.debug('KEYS:' + str(f.keys())) - raise Exception('Key not found in hdf5n file') + raise Exception('Key not found in hdf5 file') elif extension == '.mmap': @@ -1626,7 +1609,7 @@ def from_zip_file_to_movie(zipfile_name, start_end = None): start and end frame to extract @return: ''' - mov = [] + mov:List = [] print('unzipping file into movie object') if start_end is not None: num_frames = start_end[1] - start_end[0] diff --git a/caiman/base/rois.py b/caiman/base/rois.py index 8b987707a..c1bb379f8 100644 --- a/caiman/base/rois.py +++ b/caiman/base/rois.py @@ -32,6 +32,7 @@ from skimage.draw import polygon import tempfile import time +from typing import List import zipfile from ..motion_correction import tile_and_correct @@ -1260,7 +1261,7 @@ def detect_duplicates_and_subsets(binary_masks, predictions=None, r_values=None, one, two = np.unravel_index(max_idx, overlap_tmp.shape) max_val = overlap_tmp[one, two] - indeces_to_keep = [] + indeces_to_keep:List = [] indeces_to_remove = [] while max_val > 0: one, two = np.unravel_index(max_idx, overlap_tmp.shape) diff --git a/caiman/base/timeseries.py b/caiman/base/timeseries.py index bcca659e1..581709547 100644 --- a/caiman/base/timeseries.py +++ b/caiman/base/timeseries.py @@ -19,6 +19,8 @@ import tifffile import warnings +from caiman.paths import memmap_frames_filename + try: cv2.setNumThreads(0) except: @@ -112,7 +114,8 @@ def __array_finalize__(self, obj): self.file_name = getattr(obj, 'file_name', None) self.meta_data = getattr(obj, 'meta_data', None) - def save(self, file_name, to32=True, order='F',imagej=False, bigtiff=True, software='CaImAn', compress=0): + def save(self, file_name, to32=True, order='F',imagej=False, bigtiff=True, + software='CaImAn', compress=0, var_name_hdf5='mov'): """ Save the timeseries in various formats @@ -126,6 +129,9 @@ def save(self, file_name, to32=True, order='F',imagej=False, bigtiff=True, softw order: 'F' or 'C' C or Fortran order + var_name_hdf5: str + Name of hdf5 file subdirectory + Raises: Exception 'Extension Unknown' @@ -191,10 +197,12 @@ def save(self, file_name, to32=True, order='F',imagej=False, bigtiff=True, softw if self.meta_data[0] is None: savemat(file_name, {'input_arr': np.rollaxis( - input_arr, axis=0, start=3), 'start_time': self.start_time, 'fr': self.fr, 'meta_data': [], 'file_name': f_name}) + input_arr, axis=0, start=3), 'start_time': self.start_time, + 'fr': self.fr, 'meta_data': [], 'file_name': f_name}) else: savemat(file_name, {'input_arr': np.rollaxis( - input_arr, axis=0, start=3), 'start_time': self.start_time, 'fr': self.fr, 'meta_data': self.meta_data, 'file_name': f_name}) + input_arr, axis=0, start=3), 'start_time': self.start_time, + 'fr': self.fr, 'meta_data': self.meta_data, 'file_name': f_name}) elif extension in ('.hdf5', '.h5'): with h5py.File(file_name, "w") as f: @@ -203,7 +211,7 @@ def save(self, file_name, to32=True, order='F',imagej=False, bigtiff=True, softw else: input_arr = np.array(self) - dset = f.create_dataset("mov", data=input_arr) + dset = f.create_dataset(var_name_hdf5, data=input_arr) dset.attrs["fr"] = self.fr dset.attrs["start_time"] = self.start_time try: @@ -228,8 +236,7 @@ def save(self, file_name, to32=True, order='F',imagej=False, bigtiff=True, softw input_arr = np.transpose(input_arr, list(range(1, len(dims) + 1)) + [0]) input_arr = np.reshape(input_arr, (np.prod(dims), T), order='F') - fname_tot = base_name + '_d1_' + str(dims[0]) + '_d2_' + str(dims[1]) + '_d3_' + str( - 1 if len(dims) == 2 else dims[2]) + '_order_' + str(order) + '_frames_' + str(T) + '_.mmap' + fname_tot = memmap_frames_filename(base_name, dims, T, order) fname_tot = os.path.join(os.path.split(file_name)[0], fname_tot) big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=(np.uint64(np.prod(dims)), np.uint64(T)), order=order) @@ -253,7 +260,6 @@ def concatenate(*args, **kwargs): """ # todo: todocument return - obj = [] frRef = None for arg in args: for m in arg: diff --git a/caiman/cluster.py b/caiman/cluster.py index 852ebe9a8..6c6a22ff2 100644 --- a/caiman/cluster.py +++ b/caiman/cluster.py @@ -411,7 +411,8 @@ def setup_cluster(backend='multiprocessing', n_processes=None, single_thread=Fal 'A cluster is already runnning. Terminate with dview.terminate() if you want to restart.') if (platform.system() == 'Darwin') and (sys.version_info > (3, 0)): try: - if 'kernel' in get_ipython().trait_names(): # If you're on OSX and you're running under Jupyter or Spyder, + if 'kernel' in get_ipython().trait_names(): # type: ignore + # If you're on OSX and you're running under Jupyter or Spyder, # which already run the code in a forkserver-friendly way, this # can eliminate some setup and make this a reasonable approach. # Otherwise, seting VECLIB_MAXIMUM_THREADS=1 or using a different diff --git a/caiman/components_evaluation.py b/caiman/components_evaluation.py index f96c4ae50..73bd8fc23 100644 --- a/caiman/components_evaluation.py +++ b/caiman/components_evaluation.py @@ -18,6 +18,7 @@ import scipy from scipy.sparse import csc_matrix from scipy.stats import norm +from typing import Any, List import warnings from caiman.paths import caiman_datadir @@ -165,7 +166,7 @@ def find_activity_intervals(C, Npeaks=5, tB=-3, tA=10, thres=0.3): # todo todocument import peakutils K, T = np.shape(C) - L = [] + L:List = [] for i in range(K): if np.sum(np.abs(np.diff(C[i, :]))) == 0: L.append([]) @@ -204,7 +205,7 @@ def classify_components_ep(Y, A, C, b, f, Athresh=0.1, Npeaks=5, tB=-3, tA=10, t LOC = find_activity_intervals(C, Npeaks=Npeaks, tB=tB, tA=tA, thres=thres) rval = np.zeros(K) - significant_samples = [] + significant_samples:List[Any] = [] for i in range(K): if i % 200 == 0: # Show status periodically logging.info('Components evaluated:' + str(i)) @@ -414,20 +415,11 @@ def evaluate_components(Y, traces, A, C, b, f, final_frate, remove_baseline=True #%% -# FIXME xrange is python2-specific -def chunker(seq, size): - for pos in xrange(0, len(seq), size): - yield seq[pos:pos + size] -#%% - def grouper(n, iterable, fillvalue=None): "grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx" args = [iter(iterable)] * n - try: # py3 - return itertools.zip_longest(*args, fillvalue=fillvalue) - except: # py2 - return itertools.izip_longest(*args, fillvalue=fillvalue) + return itertools.zip_longest(*args, fillvalue=fillvalue) #%% @@ -556,7 +548,7 @@ def select_components_from_metrics(A, dims, gSig, r_values, comp_SNR, idx_components_r = np.where(r_values >= r_values_min)[0] idx_components_raw = np.where(comp_SNR > min_SNR)[0] - idx_components = [] + idx_components:Any = [] # changes type over the function if use_cnn: # normally 1 diff --git a/caiman/mmapping.py b/caiman/mmapping.py index 476f1a739..51e7d982f 100644 --- a/caiman/mmapping.py +++ b/caiman/mmapping.py @@ -18,12 +18,15 @@ import logging import numpy as np import os +import pickle import sys import tifffile +from typing import Any, Dict, List, Optional, Tuple, Union import caiman as cm +from caiman.paths import memmap_frames_filename -def prepare_shape(mytuple): +def prepare_shape(mytuple:Tuple) -> Tuple: """ This promotes the elements inside a shape into np.uint64. It is intended to prevent overflows with some numpy operations that are sensitive to it, e.g. np.memmap """ if not isinstance(mytuple, tuple): @@ -31,7 +34,7 @@ def prepare_shape(mytuple): return tuple(map(lambda x: np.uint64(x), mytuple)) #%% -def load_memmap(filename, mode='r'): +def load_memmap(filename:str, mode:str='r') -> Union[Tuple[Any, Tuple[int,int],int], Tuple[Any, Tuple[int,int,int], int]]: """ Load a memory mapped file created by the function save_memmap Args: @@ -71,8 +74,8 @@ def load_memmap(filename, mode='r'): raise Exception('Unknown file extension (should be .mmap)') #%% -def save_memmap_each(fnames, dview=None, base_name=None, resize_fact=(1, 1, 1), remove_init=0, - idx_xy=None, xy_shifts=None, add_to_movie=0, border_to_0=0, order = 'C', slices=None): +def save_memmap_each(fnames:List[str], dview=None, base_name:str=None, resize_fact=(1, 1, 1), remove_init:int=0, + idx_xy=None, xy_shifts=None, add_to_movie:float=0, border_to_0:int=0, order:str='C', slices=None) -> List[str]: """ Create several memory mapped files using parallel processing @@ -106,6 +109,8 @@ def save_memmap_each(fnames, dview=None, base_name=None, resize_fact=(1, 1, 1), order: (undocumented) + slices: (undocumented) + Returns: fnames_tot: list paths to the created memory map files @@ -140,7 +145,7 @@ def save_memmap_each(fnames, dview=None, base_name=None, resize_fact=(1, 1, 1), #%% -def save_memmap_join(mmap_fnames, base_name=None, n_chunks=20, dview=None, add_to_mov=0): +def save_memmap_join(mmap_fnames:List[str], base_name:str=None, n_chunks:int=20, dview=None, add_to_mov=0) -> str: """ Makes a large file memmap from a number of smaller files @@ -153,6 +158,8 @@ def save_memmap_join(mmap_fnames, base_name=None, n_chunks=20, dview=None, add_t dview: cluster handle + add_to_mov: (undocumented) + """ tot_frames = 0 @@ -171,9 +178,7 @@ def save_memmap_join(mmap_fnames, base_name=None, n_chunks=20, dview=None, add_t base_name = base_name[:base_name.find( '_d1_')] + '-#-' + str(len(mmap_fnames)) - fname_tot = (base_name + '_d1_' + str(dims[0]) + '_d2_' + str(dims[1]) + '_d3_' + - str(1 if len(dims) == 2 else dims[2]) + '_order_' + str(order) + - '_frames_' + str(tot_frames) + '_.mmap') + fname_tot = memmap_frames_filename(base_name, dims, tot_frames, order) fname_tot = os.path.join(os.path.split(mmap_fnames[0])[0], fname_tot) logging.info("Memmap file for fname_tot: " + str(fname_tot)) @@ -207,7 +212,7 @@ def save_memmap_join(mmap_fnames, base_name=None, n_chunks=20, dview=None, add_t sys.stdout.flush() return fname_tot -def my_map(dv, func, args): +def my_map(dv, func, args) -> List: v = dv rc = v.client # scatter 'id', so id=0,1,2 on engines 0,1,2 @@ -216,7 +221,7 @@ def my_map(dv, func, args): amr = v.map(func, args) pending = set(amr.msg_ids) - results_all = dict() + results_all:Dict = dict() counter = 0 while pending: try: @@ -247,7 +252,7 @@ def my_map(dv, func, args): del results_all return result_ordered -def save_portion(pars): +def save_portion(pars) -> int: # todo: todocument use_mmap_save = False big_mov, d, tot_frames, fnames, idx_start, idx_end, add_to_mov = pars @@ -285,7 +290,7 @@ def save_portion(pars): #%% -def save_place_holder(pars): +def save_place_holder(pars:List) -> str: """ To use map reduce """ # todo: todocument @@ -299,9 +304,9 @@ def save_place_holder(pars): #%% -def save_memmap(filenames, base_name='Yr', resize_fact=(1, 1, 1), remove_init=0, idx_xy=None, - order='F', xy_shifts=None, is_3D=False, add_to_movie=0, border_to_0=0, dview = None, - n_chunks=100, slices=None): +def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1, 1), remove_init:int=0, idx_xy:Tuple=None, + order:str='F', xy_shifts:Optional[List]=None, is_3D:bool=False, add_to_movie:float=0, border_to_0=0, dview = None, + n_chunks:int=100, slices=None) -> str: """ Efficiently write data from a list of tif files into a memory mappable file @@ -370,6 +375,7 @@ def save_memmap(filenames, base_name='Yr', resize_fact=(1, 1, 1), remove_init=0, logging.debug('Distributing memory map over many files') # Here we make a bunch of memmap files in the right order. Same parameters + # TODO: Use separate variables to hold the list here vs the string we return fname_new = cm.save_memmap_each(filenames, base_name = base_name, order = order, @@ -386,7 +392,7 @@ def save_memmap(filenames, base_name='Yr', resize_fact=(1, 1, 1), remove_init=0, # The goal is to make a single large memmap file, which we do here if order == 'F': - raise exception('You cannot merge files in F order, they must be in C order for CaImAn') + raise Exception('You cannot merge files in F order, they must be in C order for CaImAn') fname_new = cm.save_memmap_join(fname_new, base_name=base_name, dview=dview, n_chunks=n_chunks) @@ -411,7 +417,10 @@ def save_memmap(filenames, base_name='Yr', resize_fact=(1, 1, 1), remove_init=0, Yr = Yr[remove_init:, idx_xy[0], idx_xy[1], idx_xy[2]] else: - Yr = cm.load(f, fr=1, in_memory=True) if (isinstance(f, basestring) or isinstance(f, list)) else cm.movie(f) # TODO: Rewrite more legibly + if isinstance(f, basestring) or isinstance(f, list): + Yr = cm.load(f, fr=1, in_memory=True) + else: + Yr = cm.movie(f) if xy_shifts is not None: Yr = Yr.apply_shifts(xy_shifts, interpolation='cubic', remove_blanks=False) @@ -485,7 +494,7 @@ def save_memmap(filenames, base_name='Yr', resize_fact=(1, 1, 1), remove_init=0, #%% -def parallel_dot_product(A, b, block_size=5000, dview=None, transpose=False, num_blocks_per_run=20): +def parallel_dot_product(A:np.ndarray, b, block_size:int=5000, dview=None, transpose=False, num_blocks_per_run=20) -> np.ndarray: # todo: todocument """ Chunk matrix product between matrix and column vectors @@ -566,10 +575,9 @@ def parallel_dot_product(A, b, block_size=5000, dview=None, transpose=False, num #%% -def dot_place_holder(par): +def dot_place_holder(par:List) -> Tuple: # todo: todocument - import pickle A_name, idx_to_pass, b_, transpose = par A_, _, _ = load_memmap(A_name) b_ = pickle.loads(b_).astype(np.float32) @@ -594,7 +602,7 @@ def dot_place_holder(par): #%% def save_tif_to_mmap_online(movie_iterable, save_base_name='YrOL_', order='C', - add_to_movie=0, border_to_0=0): + add_to_movie=0, border_to_0=0) -> str: # todo: todocument if isinstance(movie_iterable, basestring): # Allow specifying a filename rather than its data rep @@ -604,11 +612,9 @@ def save_tif_to_mmap_online(movie_iterable, save_base_name='YrOL_', order='C', count = 0 new_mov = [] - dims = (len(movie_iterable),) + movie_iterable[0].shape + dims = (len(movie_iterable),) + movie_iterable[0].shape # TODO: Don't pack length into dims - fname_tot = (save_base_name + '_d1_' + str(dims[1]) + '_d2_' + str(dims[2]) + '_d3_' + - str(1 if len(dims) == 3 else dims[3]) + '_order_' + str(order) + - '_frames_' + str(dims[0]) + '_.mmap') + fname_tot = (memmap_frames_filename(save_base_name, dims[1:], dims[0], order)) big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=prepare_shape((np.prod(dims[1:]), dims[0])), order=order) diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py index 446d64e89..967869905 100644 --- a/caiman/motion_correction.py +++ b/caiman/motion_correction.py @@ -57,8 +57,12 @@ import os import pylab as pl import tifffile +from typing import List, Optional import caiman as cm +import caiman.base.movies +import caiman.motion_correction +from caiman.paths import memmap_frames_filename from .mmapping import prepare_shape try: @@ -93,7 +97,7 @@ class implementing motion correction operations def __init__(self, fname, min_mov=None, dview=None, max_shifts=(6, 6), niter_rig=1, splits_rig=14, num_splits_to_process_rig=None, strides=(96, 96), overlaps=(32, 32), splits_els=14, num_splits_to_process_els=[7, None], upsample_factor_grid=4, max_deviation_rigid=3, shifts_opencv=True, nonneg_movie=True, gSig_filt=None, - use_cuda=False, border_nan=True, pw_rigid=False, num_frames_split=80): + use_cuda=False, border_nan=True, pw_rigid=False, num_frames_split=80, var_name_hdf5='mov'): """ Constructor class for motion correction operations @@ -159,6 +163,9 @@ def __init__(self, fname, min_mov=None, dview=None, max_shifts=(6, 6), niter_rig Number of frames in each batch. Used when cosntructing the options through the params object + var_name_hdf5: str, default: 'mov' + If loading from hdf5, name of the variable to load + Returns: self @@ -190,6 +197,7 @@ def __init__(self, fname, min_mov=None, dview=None, max_shifts=(6, 6), niter_rig self.use_cuda = use_cuda self.border_nan = border_nan self.pw_rigid = pw_rigid + self.var_name_hdf5 = var_name_hdf5 if self.use_cuda and not HAS_CUDA: logging.debug("pycuda is unavailable. Falling back to default FFT.") @@ -216,10 +224,12 @@ def motion_correct(self, template=None, save_movie=False): if self.min_mov is None: if self.gSig_filt is None: self.min_mov = np.array([cm.load(self.fname[0], + var_name_hdf5=self.var_name_hdf5, subindices=slice(400))]).min() else: self.min_mov = np.array([high_pass_filter_space(m_, self.gSig_filt) - for m_ in cm.load(self.fname[0], subindices=slice(400))]).min() + for m_ in cm.load(self.fname[0], var_name_hdf5=self.var_name_hdf5, + subindices=slice(400))]).min() if self.pw_rigid: self.motion_correct_pwrigid(template=template, save_movie=save_movie) @@ -258,9 +268,9 @@ def motion_correct_rigid(self, template=None, save_movie=False): logging.debug('Entering Rigid Motion Correction') logging.debug(-self.min_mov) # XXX why the minus? self.total_template_rig = template - self.templates_rig = [] - self.fname_tot_rig = [] - self.shifts_rig = [] + self.templates_rig:List = [] + self.fname_tot_rig:List = [] + self.shifts_rig:List = [] for fname_cur in self.fname: _fname_tot_rig, _total_template_rig, _templates_rig, _shifts_rig = motion_correct_batch_rigid( @@ -277,7 +287,8 @@ def motion_correct_rigid(self, template=None, save_movie=False): nonneg_movie=self.nonneg_movie, gSig_filt=self.gSig_filt, use_cuda=self.use_cuda, - border_nan=self.border_nan) + border_nan=self.border_nan, + var_name_hdf5=self.var_name_hdf5) if template is None: self.total_template_rig = _total_template_rig @@ -328,11 +339,11 @@ def motion_correct_pwrigid( else: self.total_template_els = template - self.fname_tot_els = [] - self.templates_els = [] - self.x_shifts_els = [] - self.y_shifts_els = [] - self.coord_shifts_els = [] + self.fname_tot_els:List = [] + self.templates_els:List = [] + self.x_shifts_els:List = [] + self.y_shifts_els:List = [] + self.coord_shifts_els:List = [] for name_cur in self.fname: for num_splits_to_process in self.num_splits_to_process_els: _fname_tot_els, new_template_els, _templates_els,\ @@ -342,7 +353,7 @@ def motion_correct_pwrigid( max_deviation_rigid=self.max_deviation_rigid, splits=self.splits_els, num_splits_to_process=num_splits_to_process, num_iter=num_iter, template=self.total_template_els, shifts_opencv=self.shifts_opencv, save_movie=save_movie, nonneg_movie=self.nonneg_movie, gSig_filt=self.gSig_filt, - use_cuda=self.use_cuda, border_nan=self.border_nan) + use_cuda=self.use_cuda, border_nan=self.border_nan, var_name_hdf5=self.var_name_hdf5) if show_template: pl.imshow(new_template_els) pl.pause(.5) @@ -457,12 +468,10 @@ def apply_shift_online(movie_iterable, xy_shifts, save_base_name=None, order='F' raise Exception('Number of shifts does not match movie length!') count = 0 new_mov = [] - dims = (len(movie_iterable),) + movie_iterable[0].shape + dims = (len(movie_iterable),) + movie_iterable[0].shape # TODO: Refactor so length is either tracked separately or is last part of tuple if save_base_name is not None: - fname_tot = save_base_name + '_d1_' + str(dims[1]) + '_d2_' + str(dims[2]) + '_d3_' + str( - 1 if len(dims) == 3 else dims[3]) + '_order_' + str(order) + '_frames_' + str(dims[0]) + '_.mmap' - + fname_tot = memmap_frames_filename(save_base_name, dims[1:], dims[0], order) big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=prepare_shape((np.prod(dims[1:]), dims[0])), order=order) @@ -516,7 +525,7 @@ def motion_correct_oneP_rigid( Returns: Motion correction object ''' - min_mov = np.array([cm.motion_correction.high_pass_filter_space( + min_mov = np.array([caiman.motion_correction.high_pass_filter_space( m_, gSig_filt) for m_ in cm.load(filename[0], subindices=range(400))]).min() new_templ = None @@ -655,7 +664,7 @@ def motion_correct_online(movie_iterable, add_to_movie, max_shift_w=25, max_shif else: init_mov = movie_iterable[slice(0, init_frames_template, 1)] - dims = (len(movie_iterable),) + movie_iterable[0].shape + dims = (len(movie_iterable),) + movie_iterable[0].shape # TODO: Refactor so length is either tracked separately or is last part of tuple logging.debug("dimensions:" + str(dims)) if use_median_as_template: @@ -677,13 +686,13 @@ def motion_correct_online(movie_iterable, add_to_movie, max_shift_w=25, max_shif min_mov = 0 buffer_size_frames = 100 buffer_size_template = 100 - buffer_frames = collections.deque(maxlen=buffer_size_frames) - buffer_templates = collections.deque(maxlen=buffer_size_template) + buffer_frames:collections.deque = collections.deque(maxlen=buffer_size_frames) + buffer_templates:collections.deque = collections.deque(maxlen=buffer_size_template) max_w, max_h, min_w, min_h = 0, 0, 0, 0 big_mov = None if return_mov: - mov = [] + mov:Optional[List] = [] else: mov = None @@ -697,8 +706,7 @@ def motion_correct_online(movie_iterable, add_to_movie, max_shift_w=25, max_shif dims = (dims[0], dims[1] + min_h - max_h, dims[2] + min_w - max_w) - fname_tot = save_base_name + '_d1_' + str(dims[1]) + '_d2_' + str(dims[2]) + '_d3_' + str( - 1 if len(dims) == 3 else dims[3]) + '_order_' + str(order) + '_frames_' + str(dims[0]) + '_.mmap' + fname_tot:Optional[str] = memmap_frames_filename(save_base_name, dims[1:], dims[0], order) big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=prepare_shape((np.prod(dims[1:]), dims[0])), order=order) @@ -764,11 +772,10 @@ def motion_correct_online(movie_iterable, add_to_movie, max_shift_w=25, max_shif new_img = new_img[:, :min_w] if (save_base_name is not None) and (n_iter == (n + 1)): + big_mov[:, idx_frame] = np.reshape(new_img, np.prod(dims[1:]), order='F') # type: ignore + # mypy cannot prove that big_mov is not still None - big_mov[:, idx_frame] = np.reshape( - new_img, np.prod(dims[1:]), order='F') - - if return_mov and (n_iter == (n + 1)): + if mov is not None and (n_iter == (n + 1)): mov.append(new_img) if show_movie: @@ -783,7 +790,7 @@ def motion_correct_online(movie_iterable, add_to_movie, max_shift_w=25, max_shif if save_base_name is not None: logging.debug('Flushing memory') - big_mov.flush() + big_mov.flush() # type: ignore # mypy cannot prove big_mov is not still None del big_mov gc.collect() @@ -2145,7 +2152,7 @@ def compute_metrics_motion_correction(fname, final_size_x, final_size_y, swap_di def motion_correct_batch_rigid(fname, max_shifts, dview=None, splits=56, num_splits_to_process=None, num_iter=1, template=None, shifts_opencv=False, save_movie_rigid=False, add_to_movie=None, nonneg_movie=False, gSig_filt=None, subidx=slice(None, None, 1), use_cuda=False, - border_nan=True): + border_nan=True, var_name_hdf5='mov'): """ Function that perform memory efficient hyper parallelized rigid motion corrections while also saving a memory mappable file @@ -2199,19 +2206,19 @@ def motion_correct_batch_rigid(fname, max_shifts, dview=None, splits=56, num_spl """ corrected_slicer = slice(subidx.start, subidx.stop, subidx.step * 10) - m = cm.load(fname, subindices=corrected_slicer) + m = cm.load(fname, var_name_hdf5=var_name_hdf5, subindices=corrected_slicer) if m.shape[0] < 300: - m = cm.load(fname, subindices=corrected_slicer) + m = cm.load(fname, var_name_hdf5=var_name_hdf5, subindices=corrected_slicer) elif m.shape[0] < 500: corrected_slicer = slice(subidx.start, subidx.stop, subidx.step * 5) - m = cm.load(fname, subindices=corrected_slicer) + m = cm.load(fname, var_name_hdf5=var_name_hdf5, subindices=corrected_slicer) else: corrected_slicer = slice(subidx.start, subidx.stop, subidx.step * 30) - m = cm.load(fname, subindices=corrected_slicer) + m = cm.load(fname, var_name_hdf5=var_name_hdf5, subindices=corrected_slicer) if len(m.shape) < 3: - m = cm.load(fname) + m = cm.load(fname, var_name_hdf5=var_name_hdf5) m = m[corrected_slicer] logging.warning("Your original file was saved as a single page " + "file. Consider saving it in multiple smaller files" + @@ -2221,7 +2228,7 @@ def motion_correct_batch_rigid(fname, max_shifts, dview=None, splits=56, num_spl m = cm.movie( np.array([high_pass_filter_space(m_, gSig_filt) for m_ in m])) - template = cm.motion_correction.bin_median( + template = caiman.motion_correction.bin_median( m.motion_correct(max_shifts[1], max_shifts[0], template=None)[0]) new_templ = template @@ -2236,7 +2243,7 @@ def motion_correct_batch_rigid(fname, max_shifts, dview=None, splits=56, num_spl save_movie = False fname_tot_rig = None - res_rig = [] + res_rig:List = [] for iter_ in range(num_iter): logging.debug(iter_) old_templ = new_templ.copy() @@ -2249,7 +2256,7 @@ def motion_correct_batch_rigid(fname, max_shifts, dview=None, splits=56, num_spl dview=dview, save_movie=save_movie, base_name=os.path.split( fname)[-1][:-4] + '_rig_', subidx = subidx, num_splits=num_splits_to_process, shifts_opencv=shifts_opencv, nonneg_movie=nonneg_movie, gSig_filt=gSig_filt, - use_cuda=use_cuda, border_nan=border_nan) + use_cuda=use_cuda, border_nan=border_nan, var_name_hdf5=var_name_hdf5) new_templ = np.nanmedian(np.dstack([r[-1] for r in res_rig]), -1) if gSig_filt is not None: @@ -2259,7 +2266,7 @@ def motion_correct_batch_rigid(fname, max_shifts, dview=None, splits=56, num_spl total_template = new_templ templates = [] - shifts = [] + shifts:List = [] for rr in res_rig: shift_info, idxs, tmpl = rr templates.append(tmpl) @@ -2271,7 +2278,7 @@ def motion_correct_batch_pwrigid(fname, max_shifts, strides, overlaps, add_to_mo dview=None, upsample_factor_grid=4, max_deviation_rigid=3, splits=56, num_splits_to_process=None, num_iter=1, template=None, shifts_opencv=False, save_movie=False, nonneg_movie=False, gSig_filt=None, - use_cuda=False, border_nan=True): + use_cuda=False, border_nan=True, var_name_hdf5='mov'): """ Function that perform memory efficient hyper parallelized rigid motion corrections while also saving a memory mappable file @@ -2361,7 +2368,7 @@ def motion_correct_batch_pwrigid(fname, max_shifts, strides, overlaps, add_to_mo upsample_factor_grid=upsample_factor_grid, order='F', dview=dview, save_movie=save_movie, base_name=os.path.split(fname)[-1][:-4] + '_els_', num_splits=num_splits_to_process, shifts_opencv=shifts_opencv, nonneg_movie=nonneg_movie, gSig_filt=gSig_filt, - use_cuda=use_cuda, border_nan=border_nan) + use_cuda=use_cuda, border_nan=border_nan, var_name_hdf5=var_name_hdf5) new_templ = np.nanmedian(np.dstack([r[-1] for r in res_el]), -1) if gSig_filt is not None: @@ -2407,23 +2414,26 @@ def tile_and_correct_wrapper(params): img_name, out_fname, idxs, shape_mov, template, strides, overlaps, max_shifts,\ add_to_movie, max_deviation_rigid, upsample_factor_grid, newoverlaps, newstrides, \ - shifts_opencv, nonneg_movie, gSig_filt, is_fiji, use_cuda, border_nan = params + shifts_opencv, nonneg_movie, gSig_filt, is_fiji, use_cuda, border_nan, var_name_hdf5 = params name, extension = os.path.splitext(img_name)[:2] extension = extension.lower() shift_info = [] - if extension == '.tif' or extension == '.tiff': # check if tiff file -# with tifffile.TiffFile(img_name) as tffl: -# imgs = tffl.asarray(img_name, key=idxs) - imgs = cm.load(img_name, subindices=idxs) - - elif extension == '.sbx': # check if sbx file - imgs = cm.base.movies.sbxread(img_name, idxs[0], len(idxs)) - elif extension == '.sima' or extension == '.hdf5' or extension == '.h5': - imgs = cm.load(img_name, subindices=list(idxs)) - elif extension == '.avi': - imgs = cm.load(img_name, subindices=np.array(idxs)) +# if extension == '.tif' or extension == '.tiff': # check if tiff file +## with tifffile.TiffFile(img_name) as tffl: +## imgs = tffl.asarray(img_name, key=idxs) +# imgs = cm.load(img_name, subindices=idxs) +# +# elif extension == '.sbx': # check if sbx file +# imgs = cm.base.movies.sbxread(img_name, idxs[0], len(idxs)) +# elif extension == '.sima' or extension == '.hdf5' or extension == '.h5': +# imgs = cm.load(img_name, subindices=list(idxs), +# var_name_hdf5=var_name_hdf5) +# elif extension == '.avi': +# imgs = cm.load(img_name, subindices=np.array(idxs)) + + imgs = cm.load(img_name, subindices=idxs) mc = np.zeros(imgs.shape, dtype=np.float32) for count, img in enumerate(imgs): if count % 10 == 0: @@ -2455,7 +2465,7 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 max_shifts=(12, 12), max_deviation_rigid=3, newoverlaps=None, newstrides=None, upsample_factor_grid=4, order='F', dview=None, save_movie=True, base_name=None, subidx = None, num_splits=None, shifts_opencv=False, nonneg_movie=False, gSig_filt=None, - use_cuda=False, border_nan=True): + use_cuda=False, border_nan=True, var_name_hdf5='mov'): """ """ @@ -2464,73 +2474,78 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 extension = extension.lower() is_fiji = False - if extension == '.tif' or extension == '.tiff': # check if tiff file - with tifffile.TiffFile(fname) as tf: - T = len(tf.pages) - if T == 1: # Fiji-generated TIF - is_fiji = True - try: - T, d1, d2 = tf[0].shape - except: - T, d1, d2 = tf.asarray().shape - tf.close() - else: - d1, d2 = tf.pages[0].shape - - elif extension == '.sbx': # check if sbx file - - shape = cm.base.movies.sbxshape(name) - d1 = shape[1] - d2 = shape[0] - T = shape[2] - - elif extension == '.sima': # check if sbx file - import sima - dataset = sima.ImagingDataset.load(fname) - shape = dataset.sequences[0].shape - d1 = shape[2] - d2 = shape[3] - T = shape[0] - del dataset - elif extension == '.npy': - raise Exception('Numpy not supported at the moment') - - elif extension in ('.hdf5', '.h5'): - with h5py.File(fname) as fl: - fkeys = list(fl.keys()) - if len(fkeys)==1: - fsiz = fl[fkeys[0]].shape - elif 'mov' in fkeys: - fsiz = fl['mov'].shape - elif 'imaging' in fkeys: - fsiz = fl['imaging'].shape - else: - print(fkeys) - raise Exception('Unsupported file key') - if len(fsiz) == 3: - T, d1, d2 = fsiz - elif len(fsiz) == 5: - T, _, d1, d2, _ = fsiz - else: - print(fsiz) - raise Exception('Unsupported file shape') - - elif extension == '.avi': - cap = cv2.VideoCapture(fname) - try: - T = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)) - d2 = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) - d1 = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)) - except: - logging.debug('Roll back top opencv 2') - T = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_COUNT)) - d2 = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_WIDTH)) - d1 = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT)) - cap.release() - - else: - raise Exception( - 'Unsupported file extension for parallel motion correction') +# if extension == '.tif' or extension == '.tiff': # check if tiff file +# with tifffile.TiffFile(fname) as tf: +# T = len(tf.pages) +# if T == 1: # Fiji-generated TIF +# is_fiji = True +# try: +# T, d1, d2 = tf[0].shape +# except: +# T, d1, d2 = tf.asarray().shape +# tf.close() +# else: +# d1, d2 = tf.pages[0].shape +# +# elif extension == '.sbx': # check if sbx file +# +# shape = cm.base.movies.sbxshape(name) +# d1 = shape[1] +# d2 = shape[0] +# T = shape[2] +# +# elif extension == '.sima': # check if sbx file +# import sima +# dataset = sima.ImagingDataset.load(fname) +# shape = dataset.sequences[0].shape +# d1 = shape[2] +# d2 = shape[3] +# T = shape[0] +# del dataset +# elif extension == '.npy': +# raise Exception('Numpy not supported at the moment') +# +# elif extension in ('.hdf5', '.h5'): +# with h5py.File(fname) as fl: +# fkeys = list(fl.keys()) +# if len(fkeys)==1: +# fsiz = fl[fkeys[0]].shape +# elif var_name_hdf5 in fkeys: +# fsiz = fl[var_name_hdf5].shape +# elif 'mov' in fkeys: +# fsiz = fl['mov'].shape +# elif 'imaging' in fkeys: +# fsiz = fl['imaging'].shape +# else: +# print(fkeys) +# raise Exception('Unsupported file key') +# if len(fsiz) == 3: +# T, d1, d2 = fsiz +# elif len(fsiz) == 5: +# T, _, d1, d2, _ = fsiz +# else: +# print(fsiz) +# raise Exception('Unsupported file shape') +# +# elif extension == '.avi': +# cap = cv2.VideoCapture(fname) +# try: +# T = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)) +# d2 = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) +# d1 = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)) +# except: +# logging.debug('Roll back top opencv 2') +# T = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_COUNT)) +# d2 = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_WIDTH)) +# d1 = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT)) +# cap.release() +# +# else: +# raise Exception( +# 'Unsupported file extension for parallel motion correction') + + dims, T = cm.source_extraction.cnmf.utilities.get_file_size(fname, var_name_hdf5=var_name_hdf5) + d1, d2 = dims if type(splits) is int: if subidx is None: @@ -2557,8 +2572,7 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 if save_movie: if base_name is None: base_name = os.path.split(fname)[1][:-4] - fname_tot = base_name + '_d1_' + str(dims[0]) + '_d2_' + str(dims[1]) + '_d3_' + str( - 1 if len(dims) == 2 else dims[2]) + '_order_' + str(order) + '_frames_' + str(T) + '_.mmap' + fname_tot = memmap_frames_filename(base_name, dims, T, order) fname_tot = os.path.join(os.path.split(fname)[0], fname_tot) np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=prepare_shape(shape_mov), order=order) @@ -2570,7 +2584,8 @@ def motion_correction_piecewise(fname, splits, strides, overlaps, add_to_movie=0 for idx in idxs: pars.append([fname, fname_tot, idx, shape_mov, template, strides, overlaps, max_shifts, np.array( add_to_movie, dtype=np.float32), max_deviation_rigid, upsample_factor_grid, - newoverlaps, newstrides, shifts_opencv, nonneg_movie, gSig_filt, is_fiji, use_cuda, border_nan]) + newoverlaps, newstrides, shifts_opencv, nonneg_movie, gSig_filt, is_fiji, + use_cuda, border_nan, var_name_hdf5]) if dview is not None: logging.info('** Starting parallel motion correction **') diff --git a/caiman/paths.py b/caiman/paths.py index 11cc4a240..92b2de63b 100644 --- a/caiman/paths.py +++ b/caiman/paths.py @@ -5,8 +5,11 @@ """ import os +from typing import Tuple -def caiman_datadir(): +####### +# datadir +def caiman_datadir() -> str: """ The datadir is a user-configurable place which holds a user-modifiable copy of data the Caiman libraries need to function, alongside code demos and other things. @@ -18,5 +21,23 @@ def caiman_datadir(): else: return os.path.join(os.path.expanduser("~"), "caiman_data") -def caiman_datadir_exists(): +def caiman_datadir_exists() -> bool: return os.path.isdir(caiman_datadir()) + +###### +# memmap files +# +# Right now these are usually stored in the cwd of the script, although the basename could change that +# In the future we may consistently store these somewhere under the caiman_datadir + +def memmap_frames_filename(basename:str, dims:Tuple, frames:int, order:str='F') -> str: + # Some functions calling this have the first part of *their* dims Tuple be the number of frames. + # They *must* pass a slice to this so dims is only X, Y, and optionally Z. Frames is passed separately. + dimfield_0 = dims[0] + dimfield_1 = dims[1] + if len(dims) == 3: + dimfield_2 = dims[2] + else: + dimfield_2 = 1 + return f"{basename}_d1_{dimfield_0}_d2_{dimfield_1}_d3_{dimfield_2}_order_{order}_frames_{frames}_.mmap" + diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index d74c01411..e49676c12 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -389,9 +389,9 @@ def fit(self, images, indeces=[slice(None), slice(None)]): if self.params.get('patch', 'rf') is None and (is_sliced or 'ndarray' in str(type(images))): images = images[indeces] self.dview = None - logging.warning("Parallel processing in a single patch\ - is not available for loaded in memory or sliced\ - data.") + logging.warning("Parallel processing in a single patch " + "is not available for loaded in memory or sliced" + + " data.") T = images.shape[0] self.params.set('online', {'init_batch': T}) @@ -400,7 +400,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): Y = np.transpose(images, list(range(1, len(self.dims) + 1)) + [0]) Yr = np.transpose(np.reshape(images, (T, -1), order='F')) if np.isfortran(Yr): - raise Exception('The file is in F order, it should be in C order (see save_memmap function') + raise Exception('The file is in F order, it should be in C order (see save_memmap function)') logging.info((T,) + self.dims) @@ -434,10 +434,10 @@ def fit(self, images, indeces=[slice(None), slice(None)]): logging.info('using ' + str(self.params.get('temporal', 'block_size_temp')) + ' block_size_temp') if self.params.get('patch', 'rf') is None: # no patches - logging.debug('preprocessing ...') + logging.info('preprocessing ...') Yr = self.preprocess(Yr) if self.estimates.A is None: - logging.debug('initializing ...') + logging.info('initializing ...') self.initialize(Y) if self.params.get('patch', 'only_init'): # only return values after initialization @@ -467,7 +467,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): logging.info(('Keeping ' + str(len(idx_components)) + ' and discarding ' + str(len(idx_components_bad)))) self.estimates.C = self.estimates.C[idx_components] - self.estimates.A = self.estimates.A[:, idx_components] + self.estimates.A = self.estimates.A[:, idx_components] # type: ignore # not provable that self.initialise provides a value self.estimates.YrA = self.estimates.YrA[idx_components] self.estimates.normalize_components() @@ -490,13 +490,9 @@ def fit(self, images, indeces=[slice(None), slice(None)]): if not self.skip_refinement: logging.info('refinement...') if self.params.get('merging', 'do_merge'): - logging.debug('merging components ...') - logging.debug(self.estimates.A.shape) - logging.debug(self.estimates.C.shape) + logging.info('merging components ...') self.merge_comps(Yr, mx=50, fast_merge=True) - logging.debug(self.estimates.A.shape) - logging.debug(self.estimates.C.shape) logging.info('Updating spatial ...') self.update_spatial(Yr, use_init=False) @@ -528,8 +524,8 @@ def fit(self, images, indeces=[slice(None), slice(None)]): else: # use patches if self.params.get('patch', 'stride') is None: self.params.set('patch', {'stride': np.int(self.params.get('patch', 'rf') * 2 * .1)}) - logging.debug( - ('**** Setting the stride to 10% of 2*rf automatically:' + str(self.params.get('patch', 'stride')))) + logging.info( + ('Setting the stride to 10% of 2*rf automatically:' + str(self.params.get('patch', 'stride')))) if type(images) is np.ndarray: raise Exception( @@ -545,7 +541,7 @@ def fit(self, images, indeces=[slice(None), slice(None)]): indeces=indeces) self.estimates.bl, self.estimates.c1, self.estimates.g, self.estimates.neurons_sn = None, None, None, None - print("merging") + logging.info("merging") self.estimates.merged_ROIs = [0] @@ -555,14 +551,14 @@ def fit(self, images, indeces=[slice(None), slice(None)]): while len(self.estimates.merged_ROIs) > 0: self.merge_comps(Yr, mx=np.Inf, fast_merge=True) - logging.debug("update temporal") + logging.info("update temporal") self.update_temporal(Yr, use_init=False) self.params.set('spatial', {'se': np.ones((1,) * len(self.dims), dtype=np.uint8)}) - logging.debug('update spatial ...') + logging.info('update spatial ...') self.update_spatial(Yr, use_init=False) - logging.debug("update temporal") + logging.info("update temporal") self.update_temporal(Yr, use_init=False) else: while len(self.estimates.merged_ROIs) > 0: @@ -707,8 +703,8 @@ def deconvolve(self, p=None, method_deconvolution=None, bas_nonneg=None, self.estimates.S = np.stack([results[1][i] for i in order]) self.estimates.bl = [results[3][i] for i in order] self.estimates.c1 = [results[4][i] for i in order] - self.estimates.g = [results[5][i] for i in order] - self.estimates.neuron_sn = [results[6][i] for i in order] + self.estimates.g = [results[6][i] for i in order] + self.estimates.neurons_sn = [results[5][i] for i in order] self.estimates.lam = [results[8][i] for i in order] self.estimates.YrA = F - self.estimates.C return self @@ -758,7 +754,7 @@ def HALS4traces(self, Yr, groups=None, use_groups=False, order=None, groups = list(map(list, update_order(Ab)[0])) self.estimates.groups = groups C, noisyC = HALS4activity(Yr, Ab, Cf, groups=self.estimates.groups, order=order, - **kwargs) + **kwargs) # FIXME: this function is not defined in this scope if update_bck: if bck_non_neg: self.estimates.f = C[:self.params.get('init', 'nb')] @@ -801,7 +797,7 @@ def HALS4footprints(self, Yr, update_bck=True, num_iter=2): except(): Cf = self.estimates.C Yr = Yr - self.estimates.b.dot(self.estimates.f) - Ab = HALS4shapes(Yr, Ab, Cf, iters=num_iter) + Ab = HALS4shapes(Yr, Ab, Cf, iters=num_iter) # FIXME: this function is not defined in this scope if update_bck: self.estimates.A = scipy.sparse.csc_matrix(Ab[:, self.params.get('init', 'nb'):]) self.estimates.b = Ab[:, :self.params.get('init', 'nb')] diff --git a/caiman/source_extraction/cnmf/deconvolution.py b/caiman/source_extraction/cnmf/deconvolution.py index a640870f5..ff750c54e 100644 --- a/caiman/source_extraction/cnmf/deconvolution.py +++ b/caiman/source_extraction/cnmf/deconvolution.py @@ -988,7 +988,7 @@ def estimate_time_constant(fluor, p=2, sn=None, lags=5, fudge_factor=1.): A = scipy.linalg.toeplitz(xc[lags + np.arange(lags)], xc[lags + np.arange(p)]) - sn**2 * np.eye(lags, p) - g = np.linalg.lstsq(A, xc[lags + 1:])[0] + g = np.linalg.lstsq(A, xc[lags + 1:], rcond=None)[0] gr = np.roots(np.concatenate([np.array([1]), -g.flatten()])) gr = old_div((gr + gr.conjugate()), 2.) np.random.seed(45) # We want some variability below, but it doesn't have to be random at diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index 48a419386..949cfe152 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -13,6 +13,7 @@ import logging from .utilities import detrend_df_f from .spatial import threshold_components +from .merging import merge_iteration from ...components_evaluation import ( evaluate_components_CNN, estimate_components_quality_auto, select_components_from_metrics) @@ -831,8 +832,8 @@ def evaluate_components(self, imgs, params, dview=None): thresh_cnn_lowest=opts['cnn_lowest'], r_values_lowest=opts['rval_lowest'], min_SNR_reject=opts['SNR_lowest']) - self.idx_components = idx_components - self.idx_components_bad = idx_components_bad + self.idx_components = idx_components.astype(int) + self.idx_components_bad = idx_components_bad.astype(int) self.SNR_comp = SNR_comp self.r_values = r_values self.cnn_preds = cnn_preds @@ -880,7 +881,7 @@ def filter_components(self, imgs, params, new_dict={}, dview=None): gSig scale values for CNN classifier Returns: - self: CNMF object + self: estimates object self.idx_components: np.array indeces of accepted components self.idx_components_bad: np.array @@ -916,8 +917,127 @@ def filter_components(self, imgs, params, new_dict={}, dview=None): return self + def manual_merge(self, components, params): + ''' merge a given list of components. The indeces + of components are not pythonic, i.e., they start from 1. Moreover, + the indeces refer to the absolute indeces, i.e., the indeces before + spliting the components in accepted and rejected. If you want to e.g. + merge components 1 from idx_components and 10 from idx_components_bad + you will to set + ``` + components = [[self.idx_components[0], self.idx_components_bad[9]]] + ``` + + Args: + components: list + list of components to be merged. Each element should be a + tuple, list or np.array of the components to be merged. No + duplicates are allowed. If you're merging only one pair (or + set) of components, then use a list with a single (list) + element + params: params object + + Returns: + self: estimates object + ''' + + ln = np.sum(np.array([len(comp) for comp in components])) + ids = set.union(*[set(comp) for comp in components]) + if ln != len(ids): + raise Exception('The given list contains duplicate entries') + + p = params.temporal['p'] + nbmrg = len(components) # number of merging operations + d = self.A.shape[0] + T = self.C.shape[1] + # we initialize the values + A_merged = scipy.sparse.lil_matrix((d, nbmrg)) + C_merged = np.zeros((nbmrg, T)) + S_merged = np.zeros((nbmrg, T)) + bl_merged = np.zeros((nbmrg, 1)) + c1_merged = np.zeros((nbmrg, 1)) + sn_merged = np.zeros((nbmrg, 1)) + g_merged = np.zeros((nbmrg, p)) + merged_ROIs = [] + + for i in range(nbmrg): + merged_ROI = list(set(components[i])) + logging.info('Merging components {}'.format(merged_ROI)) + merged_ROIs.append(merged_ROI) + + Acsc = self.A.tocsc()[:, merged_ROI] + Ctmp = np.array(self.C[merged_ROI]) + np.array(self.YrA[merged_ROI]) + + C_to_norm = np.sqrt(np.ravel(Acsc.power(2).sum( + axis=0)) * np.sum(Ctmp ** 2, axis=1)) + indx = np.argmax(C_to_norm) + g_idx = [merged_ROI[indx]] + fast_merge = True + bm, cm, computedA, computedC, gm, \ + sm, ss = merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, + None, g_idx, indx, params.temporal) + + A_merged[:, i] = computedA + C_merged[i, :] = computedC + S_merged[i, :] = ss[:T] + bl_merged[i] = bm + c1_merged[i] = cm + sn_merged[i] = sm + g_merged[i, :] = gm + + empty = np.ravel((C_merged.sum(1) == 0) + (A_merged.sum(0) == 0)) + nbmrg -= len(empty) + if np.any(empty): + A_merged = A_merged[:, ~empty] + C_merged = C_merged[~empty] + S_merged = S_merged[~empty] + bl_merged = bl_merged[~empty] + c1_merged = c1_merged[~empty] + sn_merged = sn_merged[~empty] + g_merged = g_merged[~empty] + + neur_id = np.unique(np.hstack(merged_ROIs)) + nr = self.C.shape[0] + good_neurons = np.setdiff1d(list(range(nr)), neur_id) + if self.idx_components is not None: + new_indeces = list(range(len(good_neurons), + len(good_neurons) + nbmrg)) + + mapping_mat = np.zeros(nr) + mapping_mat[good_neurons] = np.arange(len(good_neurons), dtype=int) + gn_ = good_neurons.tolist() + new_idx = [mapping_mat[i] for i in self.idx_components if i in gn_] + new_idx_bad = [mapping_mat[i] for i in self.idx_components_bad if i in gn_] + new_idx.sort() + new_idx_bad.sort() + self.idx_components = np.array(new_idx + new_indeces, dtype=int) + self.idx_components_bad = np.array(new_idx_bad, dtype=int) + + self.A = scipy.sparse.hstack((self.A.tocsc()[:, good_neurons], + A_merged.tocsc())) + self.C = np.vstack((self.C[good_neurons, :], C_merged)) + # we continue for the variables + if self.S is not None: + self.S = np.vstack((self.S[good_neurons, :], S_merged)) + if self.bl is not None: + self.bl = np.hstack((self.bl[good_neurons], + np.array(bl_merged).flatten())) + if self.c1 is not None: + self.c1 = np.hstack((self.c1[good_neurons], + np.array(c1_merged).flatten())) + if self.sn is not None: + self.sn = np.hstack((self.sn[good_neurons], + np.array(sn_merged).flatten())) + if self.g is not None: + self.g = np.vstack((np.vstack(self.g)[good_neurons], g_merged)) + self.nr = nr - len(neur_id) + len(C_merged) + if self.coordinates is not None: + self.coordinates = caiman.utils.visualization.get_contours(self.A,\ + self.dims, thr_method='max', thr='0.2') + def threshold_spatial_components(self, maxthr=0.25, dview=None, **kwargs): - ''' threshold spatial components. See parameters of spatial.threshold_components + ''' threshold spatial components. See parameters of + spatial.threshold_components @param medw: @param thr_method: diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index 406619481..265d6c307 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -16,6 +16,7 @@ from builtins import range import cv2 +import logging from math import sqrt import matplotlib.animation as animation import matplotlib.pyplot as plt @@ -77,7 +78,7 @@ def downscale(Y, ds, opencv=False): else: Y_ds = movie(Y).resize(fx=1. / ds[0], fy=1. / ds[1], fz=1. / ds[2], interpolation=cv2.INTER_AREA) - print('***** OPENCV!!!!') + logging.info('Downscaling using OpenCV') else: if d > 3: # raise NotImplementedError @@ -264,7 +265,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter gSiz = np.round(np.asarray(gSiz) / ssub).astype(np.int) if normalize_init is True: - print('Noise Normalization') + logging.info('Variance Normalization') if img is None: img = np.mean(Y, axis=-1) img += np.median(img) @@ -280,11 +281,11 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter if ssub != 1 or tsub != 1: if method == 'corr_pnr': - print("Spatial Downsampling 1-photon") + logging.info("Spatial downsampling 1-photon") # this icrements the performance against ground truth and solves border problems Y_ds = downscale(Y, tuple([ssub] * len(d) + [tsub]), opencv=False) else: - print("Spatial Downsampling 2-photon") + logging.info("Spatial downsampling 2-photon") # this icrements the performance against ground truth and solves border problems Y_ds = downscale(Y, tuple([ssub] * len(d) + [tsub]), opencv=True) # mean_val = np.mean(Y) @@ -296,14 +297,14 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter if nb > min(np.prod(ds), Y_ds.shape[-1]): nb = -1 - print('Roi Extraction...') + logging.info('Roi Initialization...') if method == 'greedy_roi': Ain, Cin, _, b_in, f_in = greedyROI( Y_ds, nr=K, gSig=gSig, gSiz=gSiz, nIter=nIter, kernel=kernel, nb=nb, rolling_sum=rolling_sum, rolling_length=rolling_length) if use_hals: - print('(Hals) Refining Components...') + logging.info('Refining Components using HALS NMF iterations') Ain, Cin, b_in, f_in = hals( Y_ds, Ain, Cin, b_in, f_in, maxIter=maxIter) elif method == 'corr_pnr': @@ -327,7 +328,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter from SourceExtraction.CNMF4Dendrites import CNMF4Dendrites from SourceExtraction.AuxilaryFunctions import GetCentersData # Get initialization for components center - print(Y_ds.transpose([2, 0, 1]).shape) + # print(Y_ds.transpose([2, 0, 1]).shape) if options_local_NMF is None: raise Exception('You need to define arguments for local NMF') else: @@ -353,7 +354,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter else: print(method) - raise Exception("Unsupported method") + raise Exception("Unsupported initialization method") K = np.shape(Ain)[-1] @@ -516,11 +517,11 @@ def sparseNMF(Y_ds, nr, max_iter_snmf=500, alpha=10e2, sigma_smooth=(.5, .5, .5) m = scipy.ndimage.gaussian_filter(np.transpose( Y_ds, [2, 0, 1]), sigma=sigma_smooth, mode='nearest', truncate=truncate) if remove_baseline: - print('REMOVING BASELINE') + logging.info('REMOVING BASELINE') bl = np.percentile(m, perc_baseline, axis=0) m1 = np.maximum(0, m - bl) else: - print('NOT REMOVING BASELINE') + logging.info('NOT REMOVING BASELINE') bl = np.zeros(m.shape[1:]) m1 = m @@ -596,7 +597,7 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, """ - print("Greedy initialization of spatial and temporal components using spatial Gaussian filtering") + logging.info("Greedy initialization of spatial and temporal components using spatial Gaussian filtering") d = np.shape(Y) Y[np.isnan(Y)] = 0 med = np.median(Y, axis=-1) @@ -610,13 +611,13 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, rho = imblur(Y, sig=gSig, siz=gSiz, nDimBlur=Y.ndim - 1, kernel=kernel) if rolling_sum: - print('USING ROLLING SUM FOR INITIALIZATION....') + logging.info('Using rolling sum for initialization (RollingGreedyROI)') rolling_filter = np.ones( (rolling_length), dtype=np.float32) / rolling_length rho_s = scipy.signal.lfilter(rolling_filter, 1., rho**2) v = np.amax(rho_s, axis=-1) else: - print('USING TOTAL SUM FOR INITIALIZATION....') + logging.info('Using total sum for initialization (GreedyROI)') v = np.sum(rho**2, axis=-1) for k in range(nr): @@ -921,7 +922,7 @@ def greedyROI_corr(Y, Y_ds, max_number=None, gSiz=None, gSig=None, center_psf=Tr raise Exception( 'Either min_corr or min_pnr are None. Both of them must be real numbers.') - print('One photon initialization..') + logging.info('One photon initialization (GreedyCorr)') o = options['temporal_params'].copy() o['s_min'] = None if o['p'] > 1: @@ -942,7 +943,7 @@ def greedyROI_corr(Y, Y_ds, max_number=None, gSiz=None, gSig=None, center_psf=Tr if ring_size_factor is not None: # background according to ringmodel - print('Compute Background') + logging.info('Computing ring model background') W, b0 = compute_W(Y_ds.reshape((-1, total_frames), order='F'), A, C, (d1, d2), ring_size_factor * gSiz, ssub=ssub_B) @@ -962,14 +963,14 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p B = compute_B(b0, W, B) # "-B" B += Y_ds.reshape((-1, total_frames), order='F') # "Y-B" - print('Update Spatial') + logging.info('Updating spatial components') A, _, C, _ = caiman.source_extraction.cnmf.spatial.update_spatial_components( B, C=C, f=np.zeros((0, total_frames), np.float32), A_in=A, sn=np.sqrt(downscale((sn**2).reshape(dims, order='F'), tuple([ssub] * len(dims))).ravel() / tsub) / ssub, b_in=np.zeros((d1 * d2, 0), np.float32), dview=None, dims=(d1, d2), **options['spatial_params']) - print('Update Temporal') + logging.info('Updating temporal components') C, A = caiman.source_extraction.cnmf.temporal.update_temporal_components( B, spr.csc_matrix(A, dtype=np.float32), np.zeros((d1 * d2, 0), np.float32), @@ -982,7 +983,7 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p if max_number is not None: max_number -= A.shape[-1] if max_number is not 0: - print('Initialization again') + logging.info('Searching for more neurons in the residual') A_R, C_R, _, _, center_R = init_neurons_corr_pnr( (B - A.dot(C)).reshape(Y_ds.shape, order='F'), max_number=max_number, gSiz=gSiz, gSig=gSig, @@ -994,13 +995,13 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p C = np.concatenate((C, C_R), 0) # 1st iteration on decimated data - print('Merge Components') + logging.info('Merging components') A, C = caiman.source_extraction.cnmf.merging.merge_components( B, A, [], C, [], C, [], o, options['spatial_params'], dview=None, thr=options['merging']['merge_thr'], mx=np.Inf, fast_merge=True)[:2] A = A.astype(np.float32) C = C.astype(np.float32) - print('Update Spatial') + logging.info('Updating spatial components') A, _, C, _ = caiman.source_extraction.cnmf.spatial.update_spatial_components( B, C=C, f=np.zeros((0, total_frames), np.float32), A_in=A, sn=np.sqrt(downscale((sn**2).reshape(dims, order='F'), @@ -1008,14 +1009,14 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p b_in=np.zeros((d1 * d2, 0), np.float32), dview=None, dims=(d1, d2), **options['spatial_params']) A = A.astype(np.float32) - print('Update Temporal') + logging.info('Updating temporal components') C, A = caiman.source_extraction.cnmf.temporal.update_temporal_components( B, spr.csc_matrix(A), np.zeros((d1 * d2, 0), np.float32), C, np.zeros((0, total_frames), np.float32), dview=None, bl=None, c1=None, sn=None, g=None, **o)[:2] - print('Compute Background Again') + logging.info('Recomputing background') # background according to ringmodel W, b0 = compute_W(Y_ds.reshape((-1, total_frames), order='F'), A.toarray(), C, (d1, d2), ring_size_factor * gSiz, ssub=ssub_B) @@ -1043,19 +1044,19 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p .reshape((np.prod(dims), K), order='F')) B += Y.reshape((-1, T), order='F') # "Y-B" - print('Merge Components') + logging.info('Merging components') A, C = caiman.source_extraction.cnmf.merging.merge_components( B, A, [], C, [], C, [], o, options['spatial_params'], dview=None, thr=options['merging']['merge_thr'], mx=np.Inf, fast_merge=True)[:2] A = A.astype(np.float32) C = C.astype(np.float32) - print('Update Spatial') + logging.info('Updating spatial components') options['spatial_params']['se'] = np.ones((1,) * len((d1, d2)), dtype=np.uint8) A, _, C, _ = caiman.source_extraction.cnmf.spatial.update_spatial_components( B, C=C, f=np.zeros((0, T), np.float32), A_in=A, sn=sn, b_in=np.zeros((np.prod(dims), 0), np.float32), dview=None, dims=dims, **options['spatial_params']) - print('Update Temporal') + logging.info('Updating temporal components') C, A, b__, f__, S, bl, c1, neurons_sn, g1, YrA, lam__ = \ caiman.source_extraction.cnmf.temporal.update_temporal_components( B, spr.csc_matrix(A, dtype=np.float32), @@ -1068,16 +1069,16 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p use_NMF = True if nb == -1: - print('Return full Background') + logging.info('Returning full background') b_in = spr.eye(len(B), dtype='float32') f_in = B elif nb > 0: - print('Estimate low rank Background') + logging.info('Estimate low rank background (rank = {0})'.format(nb)) print(nb) if use_NMF: model = NMF(n_components=nb, init='nndsvdar') b_in = model.fit_transform(np.maximum(B, 0)) - #f_in = model.components_.squeeze() + # f_in = model.components_.squeeze() f_in = np.linalg.lstsq(b_in, B)[0] else: b_in, s_in, f_in = spr.linalg.svds(B, k=nb) @@ -1086,12 +1087,12 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p b_in = np.empty((A.shape[0], 0)) f_in = np.empty((0, T)) if nb == 0: - print('Return Background as b and W') + logging.info('Returning background as b0 and W') return (A, C, center.T, b_in.astype(np.float32), f_in.astype(np.float32), (S.astype(np.float32), bl, c1, neurons_sn, g1, YrA, W, b0)) else: - print("Don't Return Background") + logging.info("Not returning background") return (A, C, center.T, b_in.astype(np.float32), f_in.astype(np.float32), None if ring_size_factor is None else (S.astype(np.float32), bl, c1, neurons_sn, g1, YrA)) @@ -1475,9 +1476,9 @@ def init_neurons_corr_pnr(data, max_number=None, gSiz=15, gSig=None, break else: if num_neurons % 100 == 1: - print(num_neurons - 1, 'neurons have been initialized') + logging.info('{0} neurons have been initialized'.format(num_neurons - 1)) - print('In total, ', num_neurons, 'neurons were initialized.') + logging.info('In total, {0} neurons were initialized.'.format(num_neurons)) # A = np.reshape(Ain[:num_neurons], (-1, d1 * d2)).transpose() A = np.reshape(Ain[:num_neurons], (-1, d1 * d2), order='F').transpose() C = Cin[:num_neurons] diff --git a/caiman/source_extraction/cnmf/map_reduce.py b/caiman/source_extraction/cnmf/map_reduce.py index c0de5f0df..719e24500 100644 --- a/caiman/source_extraction/cnmf/map_reduce.py +++ b/caiman/source_extraction/cnmf/map_reduce.py @@ -91,19 +91,20 @@ def cnmf_patches(args_in): logger = logging.getLogger(__name__) name_log = os.path.basename( file_name[:-5]) + '_LOG_ ' + str(idx_[0]) + '_' + str(idx_[-1]) - #logger = logging.getLogger(name_log) - #hdlr = logging.FileHandler('./' + name_log) - #formatter = logging.Formatter('%(asctime)s %(levelname)s %(message)s') - #hdlr.setFormatter(formatter) - #logger.addHandler(hdlr) - #logger.setLevel(logging.INFO) + # logger = logging.getLogger(name_log) + # hdlr = logging.FileHandler('./' + name_log) + # formatter = logging.Formatter('%(asctime)s %(levelname)s %(message)s') + # hdlr.setFormatter(formatter) + # logger.addHandler(hdlr) + # logger.setLevel(logging.INFO) logger.debug(name_log + 'START') logger.debug(name_log + 'Read file') Yr, dims, timesteps = load_memmap(file_name) - # slicing array (takes the min and max index in n-dimensional space and cuts the box they define) + # slicing array (takes the min and max index in n-dimensional space and + # cuts the box they define) # for 2d a rectangle/square, for 3d a rectangular cuboid/cube, etc. upper_left_corner = min(idx_) lower_right_corner = max(idx_) @@ -115,7 +116,7 @@ def cnmf_patches(args_in): images = np.reshape(Yr.T, [timesteps] + list(dims), order='F') if params.get('patch', 'in_memory'): - images = np.array(images[slices],dtype=np.float32) + images = np.array(images[slices], dtype=np.float32) else: images = images[slices] @@ -132,16 +133,18 @@ def cnmf_patches(args_in): cnm = cnm.fit(images) return [idx_, shapes, scipy.sparse.coo_matrix(cnm.estimates.A), - cnm.estimates.b, cnm.estimates.C, cnm.estimates.f, cnm.estimates.S, cnm.estimates.bl, cnm.estimates.c1, - cnm.estimates.neurons_sn, cnm.estimates.g, cnm.estimates.sn, cnm.params.to_dict(), cnm.estimates.YrA] + cnm.estimates.b, cnm.estimates.C, cnm.estimates.f, + cnm.estimates.S, cnm.estimates.bl, cnm.estimates.c1, + cnm.estimates.neurons_sn, cnm.estimates.g, cnm.estimates.sn, + cnm.params.to_dict(), cnm.estimates.YrA] else: return None +# %% -#%% -def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, - border_pix=0, low_rank_background=True, del_duplicates=False, - indeces=[slice(None)]*3): +def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, + memory_fact=1, border_pix=0, low_rank_background=True, + del_duplicates=False, indeces=[slice(None)]*3): """Function that runs CNMF in patches Either in parallel or sequentially, and return the result for each. @@ -217,8 +220,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, strides = stride params_copy = deepcopy(params) - - npx_per_proc = np.int(old_div(np.prod(rfs), memory_fact)) + npx_per_proc = np.prod(rfs) // memory_fact params_copy.set('preprocess', {'n_pixels_per_process': npx_per_proc}) params_copy.set('spatial', {'n_pixels_per_process': npx_per_proc}) params_copy.set('temporal', {'n_pixels_per_process': npx_per_proc}) @@ -235,7 +237,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, foo[id_f] = 1 patch_centers.append(scipy.ndimage.center_of_mass( foo.reshape(dims, order='F'))) - print(id_2d) + logging.info('Patch size: {0}'.format(id_2d)) st = time.time() if dview is not None: if 'multiprocessing' in str(type(dview)): @@ -248,12 +250,13 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, print('Something went wrong') raise finally: - print('You may think that it went well but reality is harsh') + logging.info('Patch processing complete') else: file_res = list(map(cnmf_patches, args_in)) - print((time.time() - st)) + logging.info('Elapsed time for processing patches: \ + {0}s'.format(str(time.time() - st).split('.')[0])) # count components count = 0 count_bgr = 0 @@ -314,7 +317,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, # instead of filling in the matrices, construct lists with their non-zero # entries and coordinates - print('Transforming patches into full matrix') + logging.info('Embedding patches results into whole FOV') for fff in file_res: if fff is not None: @@ -365,7 +368,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, else: empty += 1 - print('Skipped %d Empty Patch', empty) + logging.debug('Skipped %d empty patches', empty) if count_bgr > 0: idx_tot_B = np.concatenate(idx_tot_B) b_tot = np.concatenate(b_tot) @@ -401,7 +404,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, optional_outputs['F'] = F_tot optional_outputs['mask'] = mask - print("Generating background") + logging.info("Constructing background") Im = scipy.sparse.csr_matrix( (1. / mask, (np.arange(d), np.arange(d))), dtype=np.float32) @@ -415,9 +418,9 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, elif low_rank_background is None: b = Im.dot(B_tot) f = F_tot - print("Leaving background components intact") + logging.info("Leaving background components intact") elif low_rank_background: - print("Compressing background components with a low rank NMF") + logging.info("Compressing background components with a low rank NMF") B_tot = Im.dot(B_tot) Bm = (B_tot) #f = np.r_[np.atleast_2d(np.mean(F_tot, axis=0)), @@ -427,8 +430,6 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, _ = mdl.fit_transform(F_tot).T f = mdl.components_.squeeze() f = np.atleast_2d(f) - #import pdb - #pdb.set_trace() for _ in range(100): f /= np.sqrt((f**2).sum(1)[:, None]) try: @@ -447,7 +448,8 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, # B_tot = scipy.sparse.coo_matrix(B_tot) f *= nB[:, None] else: - print('Removing overlapping background components from different patches') + logging.info('Removing overlapping background components \ + from different patches') nA = np.ravel(np.sqrt(A_tot.power(2).sum(0))) A_tot /= nA A_tot = scipy.sparse.coo_matrix(A_tot) @@ -476,8 +478,8 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, memory_fact=1, b = B_tot f = F_tot - print('******** USING ONE BACKGROUND PER PATCH ******') + logging.info('using one background component per patch') - print("Generating background DONE") + logging.info("Constructing background DONE") return A_tot, C_tot, YrA_tot, b, f, sn_tot, optional_outputs diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index 485858aa5..5b162cd67 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -245,7 +245,7 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, nr = nr - len(neur_id) + len(C_merged) else: - print('No neurons merged!') + logging.info('No more components merged!') merged_ROIs = [] return A, C, nr, merged_ROIs, S, bl, c1, sn, g @@ -266,7 +266,7 @@ def merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx, indx, temporal_ computedA = np.maximum( Acsc.dot(Ctmp.dot(computedC.T)) / (computedC * computedC.T), 0) else: - print('Simple Merging Take Best Neuron') + logging.info('Simple merging ny taking best neuron') computedC = Ctmp[indx] computedA = Acsc[:, indx] # then we de-normalize them using A_to_norm diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index 0ff681b34..a043a3785 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/caiman/source_extraction/cnmf/online_cnmf.py @@ -552,7 +552,7 @@ def initialize_online(self): fls = self.params.get('data', 'fnames') opts = self.params.get_group('online') Y = caiman.load(fls[0], subindices=slice(0, opts['init_batch'], - None)).astype(np.float32) + None), var_name_hdf5=self.params.get('data', 'var_name_hdf5')).astype(np.float32) ds_factor = np.maximum(opts['ds_factor'], 1) if ds_factor > 1: @@ -621,7 +621,7 @@ def initialize_online(self): elif self.params.get('online', 'init_method') == 'seeded': self.estimates.A, self.estimates.b, self.estimates.C, self.estimates.f, self.estimates.YrA = seeded_initialization( - Y.transpose(1, 2, 0), self.estimates.A, gnb=self.params.get('init', 'nb'), k=self.params.get('init', 'k'), + Y.transpose(1, 2, 0), self.estimates.A, gnb=self.params.get('init', 'nb'), k=self.params.get('init', 'K'), gSig=self.params.get('init', 'gSig'), return_object=False) self.estimates.S = np.zeros_like(self.estimates.C) nr = self.estimates.C.shape[0] @@ -720,7 +720,8 @@ def fit_online(self, **kwargs): for file_count, ffll in enumerate(process_files): print('Now processing file ' + ffll) - Y_ = caiman.load(ffll, subindices=slice(init_batc_iter[file_count], None, None)) + Y_ = caiman.load(ffll, var_name_hdf5=self.params.get('data', 'var_name_hdf5'), + subindices=slice(init_batc_iter[file_count], None, None)) old_comps = self.N # number of existing components for frame_count, frame in enumerate(Y_): # process each file @@ -969,7 +970,7 @@ def seeded_initialization(Y, Ain, dims=None, init_batch=1000, order_init=None, g not_px = np.array(not_px).flatten() Yr = np.reshape(Y, (Ain.shape[0], Y.shape[-1]), order='F') model = NMF(n_components=gnb, init='nndsvdar', max_iter=10) - b_temp = model.fit_transform(np.maximum(Yr[not_px], 0), iter=20) + b_temp = model.fit_transform(np.maximum(Yr[not_px], 0)) f_in = model.components_.squeeze() f_in = np.atleast_2d(f_in) Y_resf = np.dot(Yr, f_in.T) @@ -1489,7 +1490,10 @@ def get_candidate_components(sv, dims, Yres_buf, min_num_trial=3, gSig=(5, 5), if na: ain /= sqrt(na) Ain.append(ain) - Y_patch.append(Yres_buf.T[indeces, :]) if compute_corr else all_indices.append(indeces) + if compute_corr: + Y_patch.append(Yres_buf.T[indeces, :]) + else: + all_indices.append(indeces) idx.append(ind) if sniper_mode: Ain_cnn.append(ain_cnn) diff --git a/caiman/source_extraction/cnmf/params.py b/caiman/source_extraction/cnmf/params.py index 13dd586ec..d7a50e58a 100644 --- a/caiman/source_extraction/cnmf/params.py +++ b/caiman/source_extraction/cnmf/params.py @@ -33,7 +33,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), sniper_mode=False, test_both=False, thresh_CNN_noisy=0.5, thresh_fitness_delta=-50, thresh_fitness_raw=None, thresh_overlap=0.5, update_freq=200, update_num_comps=True, use_dense=True, use_peak_max=True, - only_init_patch=True, params_dict={}, + only_init_patch=True, var_name_hdf5='mov', params_dict={}, ): """Class for setting the processing parameters. All parameters for CNMF, online-CNMF, quality testing, and motion correction can be set here and then used in the various processing pipeline steps. @@ -45,7 +45,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), Args: Any parameter that is not set get a default value specified by the dictionary default options - DATA PARAMETERS (CNMFParams.data) ##### + DATA PARAMETERS (CNMFParams.data) ##### fnames: list[str] list of complete paths to files that need to be processed @@ -62,6 +62,9 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), dxy: (float, float) spatial resolution of FOV in pixels per um + var_name_hdf5: str, default: 'mov' + if loading from hdf5 name of the variable to load + PATCH PARAMS (CNMFParams.patch)###### rf: int or None, default: None @@ -523,7 +526,8 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'dims': dims, 'fr': fr, 'decay_time': decay_time, - 'dxy': dxy + 'dxy': dxy, + 'var_name_hdf5': var_name_hdf5 } self.patch = { @@ -720,11 +724,11 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), self.change_params(params_dict) if self.data['dims'] is None and self.data['fnames'] is not None: - self.data['dims'] = get_file_size(self.data['fnames'])[0] + self.data['dims'] = get_file_size(self.data['fnames'], var_name_hdf5=self.data['var_name_hdf5'])[0] if self.data['fnames'] is not None: if isinstance(self.data['fnames'], str): self.data['fnames'] = [self.data['fnames']] - T = get_file_size(self.data['fnames'])[1] + T = get_file_size(self.data['fnames'], var_name_hdf5=self.data['var_name_hdf5'])[1] if len(self.data['fnames']) > 1: T = T[0] num_splits = T//max(self.motion['num_frames_split'],10) @@ -743,7 +747,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), self.init['gSiz'] = [2*gs + 1 for gs in self.init['gSig']] if gnb <= 0: - logging.warning("gnb={}, hence setting keys nb_patch and low_rank_background ".format(gnb) + + logging.warning("gnb={0}, hence setting keys nb_patch and low_rank_background ".format(gnb) + "in group patch automatically.") self.set('patch', {'nb_patch': gnb, 'low_rank_background': None}) if gnb == -1: @@ -766,18 +770,18 @@ def set(self, group, val_dict, set_if_not_exists=False, verbose=False): """ if not hasattr(self, group): - raise KeyError('No group in CNMFParams named {}'.format(group)) + raise KeyError('No group in CNMFParams named {0}'.format(group)) d = getattr(self, group) for k, v in val_dict.items(): if k not in d and not set_if_not_exists: if verbose: logging.warning( - "NOT setting value of key {} in group {}, because no prior key existed...".format(k, group)) + "NOT setting value of key {0} in group {1}, because no prior key existed...".format(k, group)) else: if np.any(d[k] != v): logging.warning( - "Changing key {} in group {} from {} to {}".format(k, group, d[k], v)) + "Changing key {0} in group {1} from {2} to {3}".format(k, group, d[k], v)) d[k] = v def get(self, group, key): @@ -791,11 +795,11 @@ def get(self, group, key): """ if not hasattr(self, group): - raise KeyError('No group in CNMFParams named {}'.format(group)) + raise KeyError('No group in CNMFParams named {0}'.format(group)) d = getattr(self, group) if key not in d: - raise KeyError('No key {} in group {}'.format(key, group)) + raise KeyError('No key {0} in group {1}'.format(key, group)) return d[key] @@ -807,7 +811,7 @@ def get_group(self, group): """ if not hasattr(self, group): - raise KeyError('No group in CNMFParams named {}'.format(group)) + raise KeyError('No group in CNMFParams named {0}'.format(group)) return getattr(self, group) @@ -848,5 +852,5 @@ def change_params(self, params_dict, verbose=False): if k in d: flag = False if flag: - logging.warning('No parameter {} found!'.format(k)) + logging.warning('No parameter {0} found!'.format(k)) return self diff --git a/caiman/source_extraction/cnmf/pre_processing.py b/caiman/source_extraction/cnmf/pre_processing.py index 643f64b94..a1eb81e1c 100644 --- a/caiman/source_extraction/cnmf/pre_processing.py +++ b/caiman/source_extraction/cnmf/pre_processing.py @@ -21,7 +21,7 @@ import scipy import shutil import tempfile - +import logging from builtins import map from builtins import range from ...mmapping import load_memmap @@ -50,8 +50,9 @@ def interpolate_missing_data(Y): Exception 'The algorithm has not been tested with missing values (NaNs). Remove NaNs and rerun the algorithm.' """ coor = [] - print('checking if missing data') + logging.info('Checking for missing data entries (NaN)') if np.any(np.isnan(Y)): + logging.info('Interpolating missing data') for idx, row in enumerate(Y): nans = np.where(np.isnan(row))[0] n_nans = np.where(~np.isnan(row))[0] @@ -183,12 +184,12 @@ def get_noise_fft(Y, noise_range=[0.25, 0.5], noise_method='logmexp', max_num_sa cv2.setNumThreads(0) except: pass - psdx = [] + psdx_list = [] for y in Y.reshape(-1, T): dft = cv2.dft(y, flags=cv2.DFT_COMPLEX_OUTPUT).squeeze()[ :len(ind)][ind] - psdx.append(np.sum(1. / T * dft * dft, 1)) - psdx = np.reshape(psdx, Y.shape[:-1] + (-1,)) + psdx_list.append(np.sum(1. / T * dft * dft, 1)) + psdx = np.reshape(psdx_list, Y.shape[:-1] + (-1,)) else: xdft = np.fft.rfft(Y, axis=-1) xdft = xdft[..., ind[:xdft.shape[-1]]] diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 6ca31d4b2..ab09f6bf7 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -9,33 +9,36 @@ # noinspection PyCompatibility from past.builtins import basestring +from past.utils import old_div from builtins import zip from builtins import map from builtins import str from builtins import range -from past.utils import old_div + +import cv2 +import logging import numpy as np +import os +import scipy from scipy.sparse import coo_matrix, csc_matrix, csr_matrix from scipy.sparse import spdiags from scipy.linalg import eig -from scipy.ndimage.morphology import generate_binary_structure, iterate_structure from scipy.ndimage import label, binary_dilation +from scipy.ndimage.filters import median_filter +from scipy.ndimage.morphology import binary_closing +from scipy.ndimage.morphology import generate_binary_structure, iterate_structure +import shutil from sklearn.decomposition import NMF -from warnings import warn -import scipy -import time import tempfile -import os -import shutil +import time +from typing import List + from ...mmapping import load_memmap, parallel_dot_product -from scipy.ndimage.filters import median_filter -from scipy.ndimage.morphology import binary_closing -import cv2 def basis_denoising(y, c, boh, sn, id2_, px): if np.size(c) > 0: - _, _, a, _, _ = lars_regression_noise(y, c, 1, sn) + _, _, a, _, _ = lars_regression_noise(y, c, 1, sn) # FIXME Undefined function else: return (None, None, None) return a, px, id2_ @@ -163,7 +166,7 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, Exception "Failed to delete: " + folder """ - print('Initializing update of Spatial Components') + #logging.info('Initializing update of Spatial Components') if expandCore is None: expandCore = iterate_structure( @@ -177,7 +180,7 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, Y, A_in, C, f, n_pixels_per_process, nb) start_time = time.time() - print('computing the distance indicators') + logging.info('Computing support of spatial components') # we compute the indicator from distance indicator ind2_, nr, C, f, b_, A_in = computing_indicator( Y, A_in, b_in, C, f, nb, method_exp, dims, min_size, max_size, dist, expandCore, dview) @@ -193,13 +196,13 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if b_in is None: b_in = b_ - print('memmaping') + logging.info('Memory mapping') # we create a memory map file if not already the case, we send Cf, a # matrix that include background components C_name, Y_name, folder = creatememmap(Y, np.vstack((C, f)), dview) # we create a pixel group array (chunks for the cnmf)for the parrallelization of the process - print('Updating Spatial Components using lasso lars') + logging.info('Updating Spatial Components using lasso lars') cct = np.diag(C.dot(C.T)) pixel_groups = [] for i in range(0, np.prod(dims) - n_pixels_per_process + 1, n_pixels_per_process): @@ -224,14 +227,14 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, px, idxs_, a = pars A_[px, idxs_] = a - print("thresholding components") + logging.info("thresholding components") A_ = threshold_components(A_, dims, dview=dview, medw=medw, thr_method=thr_method, maxthr=maxthr, nrgthr=nrgthr, extract_cc=extract_cc, se=se, ss=ss) ff = np.where(np.sum(A_, axis=0) == 0) # remove empty components if np.size(ff) > 0: ff = ff[0] - print('eliminating {} empty spatial components'.format(len(ff))) + logging.info('eliminating {0} empty spatial components'.format(len(ff))) A_ = np.delete(A_, list(ff[ff < nr]), 1) C = np.delete(C, list(ff[ff < nr]), 0) nr = nr - len(ff[ff < nr]) @@ -245,7 +248,7 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, b_in = np.delete(b_in, background_ff, 1) A_ = A_[:, :nr] A_ = coo_matrix(A_) - print("Computing residuals") + logging.info("Computing residuals") if 'memmap' in str(type(Y)): Y_resf = parallel_dot_product(Y, f.T, dview=dview, block_size=block_size_spat, num_blocks_per_run=num_blocks_per_run_spat) - \ @@ -271,10 +274,12 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, # b = np.delete(b_in, background_ff, 0) # except NameError: b = b_in - print(("--- %s seconds ---" % (time.time() - start_time))) + # print(("--- %s seconds ---" % (time.time() - start_time))) + logging.info('Updating done in ' + + '{0}s'.format(str(time.time() - start_time).split(".")[0])) try: # clean up # remove temporary file created - print("Removing tempfiles created") + logging.info("Removing created tempfiles") shutil.rmtree(folder) except: raise Exception("Failed to delete: " + folder) @@ -853,7 +858,7 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, Coor['z'] = np.kron(list(range(d3)), np.ones(d2 * d1)) if not dist == np.inf: # determine search area for each neuron cm = np.zeros((nr, len(dims))) # vector for center of mass - Vr = [] # cell(nr,1); + Vr:List = [] # cell(nr,1); dist_indicator = [] pars = [] # for each dim diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 00ec2d990..86a7610e2 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -1,14 +1,14 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -"""A set of routines for estimating the temporal components, given the spatial components and temporal components - -@author: agiovann +"""A set of routines for estimating the temporal components, given the spatial +components and temporal components """ from builtins import str from builtins import map from builtins import range +import logging from scipy.sparse import spdiags, diags, coo_matrix # ,csgraph import scipy import numpy as np @@ -197,7 +197,7 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N C = Cin.copy() nA = np.ravel(A.power(2).sum(axis=0)) - print('Generating residuals') + logging.info('Generating residuals') # dview_res = None if block_size >= 500 else dview if 'memmap' in str(type(Y)): YA = parallel_dot_product(Y, A, dview=dview, block_size=block_size_temp, @@ -210,7 +210,7 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N YrA = YA - AA.T.dot(Cin).T # creating the patch of components to be computed in parrallel parrllcomp, len_parrllcomp = update_order_greedy(AA[:nr, :][:, :nr]) - print("entering the deconvolution ") + logging.info("entering the deconvolution ") C, S, bl, YrA, c1, sn, g, lam = update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, ITER, YrA, c1, sn, g, Cin, T, nA, dview, debug, AA, kwargs) ff = np.where(np.sum(C, axis=1) == 0) # remove empty components @@ -377,8 +377,9 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, YrA -= AA[jo, :].T.dot(Ctemp - C[jo, :]).T C[jo, :] = Ctemp.copy() S[jo, :] = Stemp - print((str(np.sum(len_parrllcomp[:count + 1])) + ' out of total ' + - str(nr) + ' temporal components updated')) + logging.info("{0} ".format(np.sum(len_parrllcomp[:count + 1])) + + "out of total {0} temporal components ".format(nr) + + "updated") for ii in np.arange(nr, nr + nb): cc = np.maximum(YrA[:, ii] + Cin[ii], -np.Inf) @@ -388,8 +389,9 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, if dview is not None and not('multiprocessing' in str(type(dview))): dview.results.clear() - if scipy.linalg.norm(Cin - C, 'fro') <= 1e-3 * scipy.linalg.norm(C, 'fro'): - print("stopping: overall temporal component not changing significantly") + if scipy.linalg.norm(Cin - C, 'fro') <= 1e-3*scipy.linalg.norm(C, 'fro'): + logging.info("stopping: overall temporal component not changing" + + " significantly") break else: # we keep Cin and do the iteration once more Cin = C diff --git a/caiman/source_extraction/cnmf/utilities.py b/caiman/source_extraction/cnmf/utilities.py index 193cd465b..34b23fc72 100644 --- a/caiman/source_extraction/cnmf/utilities.py +++ b/caiman/source_extraction/cnmf/utilities.py @@ -22,6 +22,7 @@ from builtins import str from builtins import range from past.utils import old_div + import cv2 import h5py import logging @@ -33,6 +34,8 @@ import scipy.ndimage.morphology as morph from skimage.feature.peak import _get_high_intensity_peaks import tifffile +from typing import List + from .initialization import greedyROI from ...base.rois import com from ...mmapping import parallel_dot_product, load_memmap @@ -797,7 +800,7 @@ def update_order_greedy(A, flag_AA=True): Eftychios A. Pnevmatikakis, Simons Foundation, 2017 """ K = np.shape(A)[-1] - parllcomp = [] + parllcomp:List = [] for i in range(K): new_list = True for ls in parllcomp: @@ -901,7 +904,25 @@ def normalize_AC(A, C, YrA, b, f, neurons_sn): return csc_matrix(A), C, YrA, b, f, neurons_sn -def get_file_size(file_name, var_name_hdf5=None): + +def get_file_size(file_name, var_name_hdf5='mov'): + """ Computes the dimensions of a file or a list of files without loading + it/them in memory. An exception is thrown if the files have FOVs with + different sizes + Args: + file_name: str or list + locations of file(s) in memory + + var_name_hdf5: 'str' + if loading from hdf5 name of the variable to load + + Returns: + dims: list + dimensions of FOV + + T: list + number of timesteps in each file + """ if isinstance(file_name, str): if os.path.exists(file_name): _, extension = os.path.splitext(file_name)[:2] @@ -918,7 +939,7 @@ def get_file_size(file_name, var_name_hdf5=None): dims[0] = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) dims[1] = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)) except(): - print('Roll back top opencv 2') + print('Roll back to opencv 2') T = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_COUNT)) dims[0] = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_WIDTH)) dims[1] = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT)) @@ -926,12 +947,12 @@ def get_file_size(file_name, var_name_hdf5=None): filename = os.path.split(file_name)[-1] Yr, dims, T = load_memmap(os.path.join( os.path.split(file_name)[0], filename)) - elif extension == '.h5' or extension == '.hdf5': + elif extension in ('.h5', '.hdf5', '.nwb'): with h5py.File(file_name, "r") as f: kk = list(f.keys()) if len(kk) == 1: siz = f[kk[0]].shape - elif var_name_hdf5 in kk: + elif var_name_hdf5 in f: siz = f[var_name_hdf5].shape else: print(kk) diff --git a/caiman/tests/test_general.py b/caiman/tests/test_general.py index dad47694d..d8a5c3663 100644 --- a/caiman/tests/test_general.py +++ b/caiman/tests/test_general.py @@ -317,8 +317,8 @@ def test_general(): logging.error("you need to set the same movie parameters than the ground truth to have a real comparison (use the comp.see() function to explore it)") pb = True if (comp.information['differences']['params_cnm']): - logging.error("you need to set the same cnmf parameters than the ground truth to have a real comparison (use the comp.see() function to explore it)") - pb = True + logging.warning("you need to set the same cnmf parameters than the ground truth to have a real comparison (use the comp.see() function to explore it)") + # pb = True if (comp.information['diff']['rig']['isdifferent']): logging.error("the rigid shifts are different from the groundtruth ") pb = True diff --git a/caiman/utils/utils.py b/caiman/utils/utils.py index d283a4446..efb4b4264 100644 --- a/caiman/utils/utils.py +++ b/caiman/utils/utils.py @@ -21,15 +21,16 @@ #\pre none import cv2 -import logging import h5py - +import logging +import matplotlib.pyplot as plt import numpy as np import os import pickle import scipy from scipy.ndimage.filters import gaussian_filter from tifffile import TiffFile +from typing import Any, Dict, List, Tuple, Union try: cv2.setNumThreads(0) @@ -37,15 +38,15 @@ pass from urllib.request import urlopen + from ..external.cell_magic_wand import cell_magic_wand from ..source_extraction.cnmf.spatial import threshold_components - from caiman.paths import caiman_datadir #%% -def download_demo(name='Sue_2x_3000_40_-46.tif', save_folder=''): +def download_demo(name:str='Sue_2x_3000_40_-46.tif', save_folder:str='') -> str: """download a file from the file list with the url of its location @@ -57,7 +58,8 @@ def download_demo(name='Sue_2x_3000_40_-46.tif', save_folder=''): save_folder: str folder inside ./example_movies to which the files will be saved. Will be created if it doesn't exist - + Returns: + Path of the saved file Raise: WrongFolder Exception """ @@ -126,8 +128,8 @@ def val_parse(v): return v -def si_parse(imd): - """parse image_description field embedded by scanimage from get iamge description +def si_parse(imd:str) -> Dict: + """parse image_description field embedded by scanimage from get image description Args: imd: image description @@ -137,15 +139,15 @@ def si_parse(imd): """ - imd = imd.split('\n') - imd = [i for i in imd if '=' in i] - imd = [i.split('=') for i in imd] - imd = [[ii.strip(' \r') for ii in i] for i in imd] - imd = {i[0]: val_parse(i[1]) for i in imd} - return imd + imddata:Any = imd.split('\n') + imddata = [i for i in imddata if '=' in i] + imddata = [i.split('=') for i in imddata] + imddata = [[ii.strip(' \r') for ii in i] for i in imddata] + imddata = {i[0]: val_parse(i[1]) for i in imddata} + return imddata -def get_image_description_SI(fname): +def get_image_description_SI(fname:str) -> List: """Given a tif file acquired with Scanimage it returns a dictionary containing the information in the image description field Args: @@ -160,8 +162,7 @@ def get_image_description_SI(fname): for idx, pag in enumerate(tf.pages): if idx % 1000 == 0: - logging.debug(idx) - # i2cd=si_parse(pag.tags['image_description'].value)['I2CData'] + logging.debug(idx) # progress report to the user field = pag.tags['image_description'].value image_descriptions.append(si_parse(field)) @@ -170,9 +171,9 @@ def get_image_description_SI(fname): #%% Generate data -def gen_data(dims=(48, 48), N=10, sig=(3, 3), tau=1., noise=.3, T=2000, - framerate=30, firerate=.5, seed=3, cmap=False, truncate=np.exp(-2), - difference_of_Gaussians=True, fluctuating_bkgrd=[50, 300]): +def gen_data(dims:Tuple[int,int]=(48, 48), N:int=10, sig:Tuple[int,int]=(3, 3), tau:float=1., noise:float=.3, T:int=2000, + framerate:int=30, firerate:float=.5, seed:int=3, cmap:bool=False, truncate:float=np.exp(-2), + difference_of_Gaussians:bool=True, fluctuating_bkgrd:List=[50, 300]) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, Tuple[int, int]]: bkgrd = 10 # fluorescence baseline np.random.seed(seed) boundary = 4 @@ -238,7 +239,6 @@ def gen_data(dims=(48, 48), N=10, sig=(3, 3), tau=1., noise=.3, T=2000, * (np.prod(dims), T)).astype('float32') + trueA.dot(trueC) if cmap: - import matplotlib.pyplot as plt import caiman as cm Y = np.reshape(Yr, dims + (T,), order='F') Cn = cm.local_correlations(Y) @@ -264,23 +264,23 @@ def gen_data(dims=(48, 48), N=10, sig=(3, 3), tau=1., noise=.3, T=2000, plt.imshow(Cn, cmap=cmap) plt.title('Correlation') plt.show() - return Yr, trueC, trueS, trueA, trueb, truef, centers, dims + return Yr, trueC, trueS, trueA, trueb, truef, centers, dims # XXX dims is always the same as passed into the function? #%% -def save_object(obj, filename): +def save_object(obj, filename:str) -> None: with open(filename, 'wb') as output: pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL) -def load_object(filename): +def load_object(filename:str) -> Any: with open(filename, 'rb') as input_obj: obj = pickle.load(input_obj) return obj #%% def apply_magic_wand(A, gSig, dims, A_thr=None, coms=None, dview=None, min_frac=0.7, max_frac=1.0, roughness=2, zoom_factor=1, - center_range=2): + center_range=2) -> np.ndarray: """ Apply cell magic Wand to results of CNMF to ease matching with labels Args: @@ -351,7 +351,7 @@ def cell_magic_wand_wrapper(params): #%% From https://codereview.stackexchange.com/questions/120802/recursively-save-python-dictionaries-to-hdf5-files-using-h5py -def save_dict_to_hdf5(dic, filename): +def save_dict_to_hdf5(dic:Dict, filename:str) -> None: ''' Save dictionary to hdf5 file Args: dic: dictionary @@ -363,7 +363,7 @@ def save_dict_to_hdf5(dic, filename): with h5py.File(filename, 'w') as h5file: recursively_save_dict_contents_to_group(h5file, '/', dic) -def load_dict_from_hdf5(filename): +def load_dict_from_hdf5(filename:str) -> Dict: ''' Load dictionary from hdf5 file Args: @@ -377,7 +377,7 @@ def load_dict_from_hdf5(filename): return recursively_load_dict_contents_from_group(h5file, '/') -def recursively_save_dict_contents_to_group(h5file, path, dic): +def recursively_save_dict_contents_to_group(h5file:h5py.File, path:str, dic:Dict) -> None: ''' Args: h5file: hdf5 object @@ -450,7 +450,7 @@ def recursively_save_dict_contents_to_group(h5file, path, dic): raise ValueError("Cannot save %s type for key '%s'." % (type(item), key)) -def recursively_load_dict_contents_from_group( h5file, path): +def recursively_load_dict_contents_from_group(h5file:h5py.File, path:str) -> Dict: '''load dictionary from hdf5 object Args: h5file: hdf5 object @@ -459,7 +459,7 @@ def recursively_load_dict_contents_from_group( h5file, path): path within the hdf5 file ''' - ans = {} + ans:Dict = {} for key, item in h5file[path].items(): if isinstance(item, h5py._hl.dataset.Dataset): diff --git a/caiman/utils/visualization.py b/caiman/utils/visualization.py index 9ab50af36..12609ec91 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -28,6 +28,7 @@ from skimage.measure import find_contours import sys from tempfile import NamedTemporaryFile +from typing import Dict from warnings import warn from ..base.rois import com @@ -230,7 +231,9 @@ def nb_view_patches(Yr, A, C, b, f, d1, d2, YrA=None, image_neurons=None, thr=0. title="Neuron Number", callback=callback) xr = Range1d(start=0, end=image_neurons.shape[1]) yr = Range1d(start=image_neurons.shape[0], end=0) - plot1 = bpl.figure(x_range=xr, y_range=yr, plot_width=300, plot_height=300) + plot1 = bpl.figure(x_range=xr, y_range=yr, + plot_width=int(min(1, d2/d1))*300, + plot_height=int(min(1, d1/d2))*300) plot1.image(image=[image_neurons[::-1, :]], x=0, y=image_neurons.shape[0], dw=d2, dh=d1, palette=grayp) @@ -286,7 +289,7 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False): # for each patches for i in range(nr): - pars = dict() + pars:Dict = dict() # we compute the cumulative sum of the energy of the Ath component that has been ordered from least to highest patch_data = A.data[A.indptr[i]:A.indptr[i + 1]] indx = np.argsort(patch_data)[::-1] diff --git a/caimanmanager.py b/caimanmanager.py index 4d1cea3c6..bf9926dbf 100755 --- a/caimanmanager.py +++ b/caimanmanager.py @@ -3,11 +3,13 @@ import argparse import distutils.dir_util import filecmp +import glob import os import shutil import string import subprocess import sys # for sys.prefix +from typing import Dict, List, Tuple from caiman.paths import caiman_datadir @@ -38,7 +40,7 @@ ############### # commands -def do_install_to(targdir, inplace=False, force=False): +def do_install_to(targdir:str, inplace:bool=False, force:bool=False) -> None: global sourcedir_base if os.path.isdir(targdir) and not force: raise Exception(targdir + " already exists") @@ -60,7 +62,7 @@ def do_install_to(targdir, inplace=False, force=False): shutil.copy(extrafile, targdir) print("Installed " + targdir) -def do_check_install(targdir, inplace=False): +def do_check_install(targdir:str, inplace:bool=False) -> None: global sourcedir_base if inplace: sourcedir_base = os.getcwd() @@ -81,8 +83,10 @@ def do_check_install(targdir, inplace=False): ok = False if ok: print("OK") + if not ok: + raise Exception("Install is dirty") -def do_run_nosetests(targdir): +def do_run_nosetests(targdir:str) -> None: out, err, ret = runcmd(["nosetests", "--traverse-namespace", "caiman"]) if ret != 0: print("Nosetests failed with return code " + str(ret)) @@ -90,7 +94,7 @@ def do_run_nosetests(targdir): else: print("Nosetests success!") -def do_run_demotests(targdir): +def do_run_demotests(targdir:str) -> None: out, err, ret = runcmd([os.path.join(caiman_datadir(), "test_demos.sh")]) if ret != 0: print("Demos failed with return code " + str(ret)) @@ -98,7 +102,7 @@ def do_run_demotests(targdir): else: print("Demos success!") -def do_nt_run_demotests(targdir): +def do_nt_run_demotests(targdir:str) -> None: # Windows platform can't run shell scripts, and doing it in batch files # is a terrible idea. So we'll do a minimal implementation of run_demos for # windows inline here. @@ -121,7 +125,7 @@ def do_nt_run_demotests(targdir): ############### # -def comparitor_all_diff_files(comparitor, path_prepend): +def comparitor_all_diff_files(comparitor, path_prepend:str): ret = list(map(lambda x: os.path.join(path_prepend, x), comparitor.diff_files)) # Initial for dirname in comparitor.subdirs.keys(): to_append = comparitor_all_diff_files(comparitor.subdirs[dirname], os.path.join(path_prepend, dirname)) @@ -129,7 +133,7 @@ def comparitor_all_diff_files(comparitor, path_prepend): ret.append(*to_append) return ret -def comparitor_all_left_only_files(comparitor, path_prepend): +def comparitor_all_left_only_files(comparitor, path_prepend:str): ret = list(map(lambda x: os.path.join(path_prepend, x), comparitor.left_only)) # Initial for dirname in comparitor.subdirs.keys(): to_append = comparitor_all_left_only_files(comparitor.subdirs[dirname], os.path.join(path_prepend, dirname)) @@ -139,7 +143,7 @@ def comparitor_all_left_only_files(comparitor, path_prepend): ############### -def runcmd(cmdlist, ignore_error=False, verbose=True): +def runcmd(cmdlist:List[str], ignore_error:bool=False, verbose:bool=True) -> Tuple[str, str, int]: # In most of my codebases, runcmd saves and returns the output. # Here I've modified it to send right to stdout, because nothing # uses the output and because the demos sometimes have issues diff --git a/demos/general/demo_pipeline.py b/demos/general/demo_pipeline.py index 704677285..7f2f36351 100755 --- a/demos/general/demo_pipeline.py +++ b/demos/general/demo_pipeline.py @@ -154,7 +154,6 @@ def main(): c, dview, n_processes = cm.cluster.setup_cluster( backend='local', 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 diff --git a/demos/general/demo_pipeline_cnmfE.py b/demos/general/demo_pipeline_cnmfE.py index ceadf0ea2..02452d86b 100755 --- a/demos/general/demo_pipeline_cnmfE.py +++ b/demos/general/demo_pipeline_cnmfE.py @@ -183,7 +183,11 @@ def main(): # %% 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, gSig=gSig[0], swap_dim=False) + cn_filter, pnr = cm.summary_images.correlation_pnr(images[::1], gSig=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 diff --git a/demos/notebooks/demo_pipeline_cnmfE.ipynb b/demos/notebooks/demo_pipeline_cnmfE.ipynb index aadded015..0313fd0c0 100644 --- a/demos/notebooks/demo_pipeline_cnmfE.ipynb +++ b/demos/notebooks/demo_pipeline_cnmfE.ipynb @@ -17,7 +17,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "#!/usr/bin/env python\n", @@ -66,7 +68,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "fnames = ['data_endoscope.tif'] # filename to be processed\n", @@ -84,7 +88,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "#%% start a cluster for parallel processing (if a cluster already exists it will be closed and a new session will be opened)\n", @@ -106,6 +112,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "collapsed": true, "scrolled": true }, "outputs": [], @@ -156,7 +163,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "if motion_correct:\n", @@ -193,7 +202,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# load memory mappable file\n", @@ -213,7 +224,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# parameters for source extraction and deconvolution\n", @@ -276,17 +289,22 @@ "metadata": {}, "source": [ "### Inspect summary images and set parameters\n", - "Check the optimal values of `min_corr` and `min_pnr` by moving slider in the figure that pops up. You can modify them in the `params` object." + "Check the optimal values of `min_corr` and `min_pnr` by moving slider in the figure that pops up. You can modify them in the `params` object. \n", + "Note that computing the correlation pnr image can be computationally and memory demanding for large datasets. In this case you can compute\n", + "only on a subset of the data (the results will not change). You can do that by changing `images[::1]` to `images[::5]` or something similar.\n", + "This will compute the correlation pnr image" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# compute some summary images (correlation and peak to noise)\n", - "cn_filter, pnr = cm.summary_images.correlation_pnr(images, gSig=gSig[0], swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile\n", + "cn_filter, pnr = cm.summary_images.correlation_pnr(images[::1], gSig=gSig[0], swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile\n", "# inspect the summary images and set the parameters\n", "inspect_correlation_pnr(cn_filter, pnr)" ] @@ -294,7 +312,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# print parameters set above, modify them if necessary based on summary images\n", @@ -312,7 +332,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "cnm = cnmf.CNMF(n_processes=n_processes, dview=dview, Ain=Ain, params=opts)\n", @@ -330,7 +352,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# cnm1 = cnmf.CNMF(n_processes, params=opts, dview=dview)\n", @@ -356,7 +380,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "#%% COMPONENT EVALUATION\n", @@ -388,7 +414,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "#%% plot contour plots of accepted and rejected components\n", @@ -406,7 +434,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# accepted components\n", @@ -417,7 +447,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# rejected components\n", @@ -435,7 +467,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "cm.stop_server(dview=dview)" @@ -452,7 +486,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# with background \n", @@ -463,7 +499,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "# without background\n", @@ -474,7 +512,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [] } @@ -482,7 +522,7 @@ "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python [default]", "language": "python", "name": "python3" }, @@ -496,7 +536,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.5.2" } }, "nbformat": 4, diff --git a/docs/Makefile b/docs/Makefile index 0ea9b027b..b7b6ef232 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -3,7 +3,8 @@ # You can set these variables from the command line. SPHINXOPTS = -SPHINXBUILD = sphinx-build -v +SPHINXBUILD = sphinx-build +SPHINXOPTS = -v PAPER = BUILDDIR = build @@ -52,38 +53,38 @@ clean: rm -rf $(BUILDDIR)/* html: - $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html + $(SPHINXBUILD) $(SPHINXOPTS) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html @echo @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." dirhtml: - $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml + $(SPHINXBUILD) $(SPHINXOPTS) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml @echo @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml." singlehtml: - $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml + $(SPHINXBUILD) $(SPHINXOPTS) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml @echo @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml." pickle: - $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle + $(SPHINXBUILD) $(SPHINXOPTS) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle @echo @echo "Build finished; now you can process the pickle files." json: - $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json + $(SPHINXBUILD) $(SPHINXOPTS) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json @echo @echo "Build finished; now you can process the JSON files." htmlhelp: - $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp + $(SPHINXBUILD) $(SPHINXOPTS) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp @echo @echo "Build finished; now you can run HTML Help Workshop with the" \ ".hhp project file in $(BUILDDIR)/htmlhelp." qthelp: - $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp + $(SPHINXBUILD) $(SPHINXOPTS) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp @echo @echo "Build finished; now you can run "qcollectiongenerator" with the" \ ".qhcp project file in $(BUILDDIR)/qthelp, like this:" @@ -92,7 +93,7 @@ qthelp: @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/Casourceseparation.qhc" applehelp: - $(SPHINXBUILD) -b applehelp $(ALLSPHINXOPTS) $(BUILDDIR)/applehelp + $(SPHINXBUILD) $(SPHINXOPTS) -b applehelp $(ALLSPHINXOPTS) $(BUILDDIR)/applehelp @echo @echo "Build finished. The help book is in $(BUILDDIR)/applehelp." @echo "N.B. You won't be able to view it unless you put it in" \ @@ -100,7 +101,7 @@ applehelp: "bundle." devhelp: - $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp + $(SPHINXBUILD) $(SPHINXOPTS) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp @echo @echo "Build finished." @echo "To view the help file:" @@ -109,84 +110,84 @@ devhelp: @echo "# devhelp" epub: - $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub + $(SPHINXBUILD) $(SPHINXOPTS) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub @echo @echo "Build finished. The epub file is in $(BUILDDIR)/epub." latex: - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + $(SPHINXBUILD) $(SPHINXOPTS) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." @echo "Run \`make' in that directory to run these through (pdf)latex" \ "(use \`make latexpdf' here to do that automatically)." latexpdf: - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + $(SPHINXBUILD) $(SPHINXOPTS) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo "Running LaTeX files through pdflatex..." $(MAKE) -C $(BUILDDIR)/latex all-pdf @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." latexpdfja: - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + $(SPHINXBUILD) $(SPHINXOPTS) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo "Running LaTeX files through platex and dvipdfmx..." $(MAKE) -C $(BUILDDIR)/latex all-pdf-ja @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." text: - $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text + $(SPHINXBUILD) $(SPHINXOPTS) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text @echo @echo "Build finished. The text files are in $(BUILDDIR)/text." man: - $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man + $(SPHINXBUILD) $(SPHINXOPTS) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man @echo @echo "Build finished. The manual pages are in $(BUILDDIR)/man." texinfo: - $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + $(SPHINXBUILD) $(SPHINXOPTS) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo @echo @echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo." @echo "Run \`make' in that directory to run these through makeinfo" \ "(use \`make info' here to do that automatically)." info: - $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + $(SPHINXBUILD) $(SPHINXOPTS) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo @echo "Running Texinfo files through makeinfo..." make -C $(BUILDDIR)/texinfo info @echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo." gettext: - $(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale + $(SPHINXBUILD) $(SPHINXOPTS) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale @echo @echo "Build finished. The message catalogs are in $(BUILDDIR)/locale." changes: - $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes + $(SPHINXBUILD) $(SPHINXOPTS) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes @echo @echo "The overview file is in $(BUILDDIR)/changes." linkcheck: - $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck + $(SPHINXBUILD) $(SPHINXOPTS) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck @echo @echo "Link check complete; look for any errors in the above output " \ "or in $(BUILDDIR)/linkcheck/output.txt." doctest: - $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest + $(SPHINXBUILD) $(SPHINXOPTS) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest @echo "Testing of doctests in the sources finished, look at the " \ "results in $(BUILDDIR)/doctest/output.txt." coverage: - $(SPHINXBUILD) -b coverage $(ALLSPHINXOPTS) $(BUILDDIR)/coverage + $(SPHINXBUILD) $(SPHINXOPTS) -b coverage $(ALLSPHINXOPTS) $(BUILDDIR)/coverage @echo "Testing of coverage in the sources finished, look at the " \ "results in $(BUILDDIR)/coverage/python.txt." xml: - $(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml + $(SPHINXBUILD) $(SPHINXOPTS) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml @echo @echo "Build finished. The XML files are in $(BUILDDIR)/xml." pseudoxml: - $(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml + $(SPHINXBUILD) $(SPHINXOPTS) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml @echo @echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml." diff --git a/docs/footer.html b/docs/footer.html index bcd108b73..bac31155f 100644 --- a/docs/footer.html +++ b/docs/footer.html @@ -1,3 +1,3 @@ \ No newline at end of file +© Copyright 2019, Flatiron Institute, Simons Foundation, New York, NY. Created using +Sphinx diff --git a/docs/source/CaImAn-Tips.md.rst b/docs/source/CaImAn-Tips.md.rst new file mode 100644 index 000000000..91f351db9 --- /dev/null +++ b/docs/source/CaImAn-Tips.md.rst @@ -0,0 +1,145 @@ +Motion Correction tips +---------------------- + +- Non-rigid motion correction is not always necessary. Sometimes, rigid + motion correction will be sufficient and it will lead to significant + performance gains in terms of speed. Check your data to before/after + rigid motion correction to decide what is best for you. The boolean + parameter ``params.motion.pw_rigid`` can be used to alternate between + rigid and non-rigid motion correction using the same function + ``MotionCorrect.motion_correct``. + +- When using piecewise rigid motion correction, use parameters that + physically make sense. For example, a typical patch size could be + around 100um x 100um since motion can in many times be approximated + as rigid for smaller patches (if the imaging is not too slow). + Similarly, the maximum allowed shifts can in typical 2p recordings + chosen to correspond to 10um. The patch size is given by the sum of + the parameters ``params.motion.strides + params.motion.overlaps``. + The maximum shifts parameter is ``params.motion.max_shifts``. These + values corresponds to pixels so make sure you have a rough idea of + the spatial resolution of your data. There is a parameter for that + ``params.data.dxy``. + +- Motion correction works in parallel by splitting each file in + multiple chunks and processing them in parallel. Make sure that the + length of each chunk is not too small by setting the parameter + ``params.motion.num_frames_split``. + +CaImAn Online Processing tips +----------------------------- + +- Important parameters for online processing are the CNN threshold + value ``params.online.thresh_CNN_noisy``, the trace SNR + ``params.online.min_SNR`` and the number of candidate components to + be considered at each timestep ``params.online.min_num_trial``. Lower + values for the thresholds (e.g., 1 for ``params.online.min_SNR`` and + 0.5 for ``params.online.thresh_CNN_noisy``) and/or higher values for + ``params.online.min_num_trial`` (e.g., 10) can lead to higher recall + values, although potentially at the expense of lower precision. In + general they are preferable for datasets that are relatively short + (e.g., 10000 frames or less). On the other hand, higher threshold + values (e.g., 1.5 for ``params.online.min_SNR`` and 0.7 for + ``params.online.thresh_CNN_noisy``) and/or lower values for + ``params.online.min_num_trial`` (e.g., 5) will lead to higher + precision values, although potentially at the expense of lower + recall. n general they are preferable for datasets that are longer + (e.g., 10000 frames or more). + +- If your analysis setup allows it, multiple epochs over the data can + be very beneficial, especially in the strict regime or high + acceptance thresholds. + +- In general, ``bare`` initialization can be used most of the times, to + capture the neuropil activity and a small number of neurons at an + initial chunk. For a large FOV with lots of active neurons, e.g., a + plane from a zebrafish dataset, ``bare`` initialization can be + inadequate. In this case, a proper initialization with ``cnmf`` can + lead to substantially better results. + +- Spatial downsampling can lead to significant speed gains, often at no + expense in terms of accuracy. It can be set through the parameter + ``ds_factor``. + +- When using the CNN for screening candidate components, the usage of a + GPU can lead to significant computational gains. + +CaImAn Batch processing tips +---------------------------- + +- In order to optimize memory consumption and parallelize computing, it + is suggested to adopt computing in patches (see companion paper). The + user will inspect the correlation image and select an appropriate + number of neurons per each patch. The + ``params.patches['rf']' and``\ params.patches.stride’ parameters + controls the size of patches and their overlap. Given the patch size + and the correlation image the user can set an upper bound on the + number of neurons per patches. We suggest to start exploring regions + that contain 4-10 neurons. + +- Important parameters for selecting components based on quality are + +- the CNN lower bound and upper threshold ``params.quality.cnn_lowest`` + and ``params.quality.min_cnn_thr`` + +- the trace SNR ``params.quality.min_SNR`` + +- the footprint consistency threshold ``params.quality.rval_thr`` + +The user should explore these parameters around the default to optimize +for specific data sets. + +1p processing tips +------------------ + +- For microendoscopic 1p data use CNMF-E’s background model and + initialization method by setting ``center_psf=True``, + ``method_init='corr_pnr'`` and ``ring_size_factor`` to some value + around 1.5. In this case the spatial and temporal components are + updated during the initialization phase, hence use + ``only_init_patch=True``. + +- Other important parameters for microendoscopic 1p data are ``gSig``, + ``gSiz``, ``min_corr`` and ``min_pnr``. ``gSig`` specifies the + gaussian width of a 2D gaussian kernel, which approximates a neuron + and ``gSiz`` the average diameter of a neuron, in general + ``4*gSig+1``. To pick the thresholds ``min_corr`` and ``min_pnr`` you + can use ``caiman.utils.visualization.inspect_correlation_pnr`` and + vary the slider values. + +- Because the background has no high spatial frequency components, it + can be spatially downscaled to speed up processing without loss in + accuracy, e.g. by a factor of 2 by setting ``ssub_B=2``. + +- The exact background can be returned as full rank matrix + (``gnb=-1``), or more compactly as parameters of the ring model + (``gnb=0``), or not at all (``gnb<-1``). Further the background can + also be approximated as low rank matrix by setting ``gnb`` to the + desired rank. ``gnb=0`` is usually the desired choice. If you have + plenty of RAM and process in patches ``gnb=-1`` is a good and faster + option. + +- The CNMF-E algorithm poses high demands on RAM. There is however a + trade off between computing time and memory usage when processing in + patches. The number of processes ``n_processes`` specifies how many + patches are processed in parallel, thus a higher number decreases + computing time but increases RAM usage. If you have insufficient RAM, + use a smaller value for ``n_processes`` to reduce memory consumption, + or don’t even use parallelization at all by setting ``dview=None``. + +Deconvolution tips +------------------ + +- Simultaneous deconvolution and source extraction can mostly offer + benefits in particularly low SNR data. In most cases, running source + extraction without deconvolution (``p=0``), followed by deconvolution + will be sufficient. + +- It is generally better to perform some sort of de-trending on the + extracted calcium traces prior to deconvolution to correct for + baseline drifts that can results in wrongfully deconvolved neural + activity. You can use the ``estimates.detrend_df_f`` methods for + that. + +- For interpreting the deconvolved neural activity varible ``S``, see + `here `__. diff --git a/docs/source/CaImAn-features-and-references.md b/docs/source/CaImAn-features-and-references.md new file mode 100644 index 000000000..696517bbc --- /dev/null +++ b/docs/source/CaImAn-features-and-references.md @@ -0,0 +1,78 @@ +## Features + +CaImAn includes a variety of scalable methods for the analysis of calcium imaging data: + +* **Handling of very large datasets** + * Memory mapping + * Parallel processing in patches + * Frame-by-frame online processing [[6]](#onacid) + * OpenCV-based efficient movie playing and resizing + +* **Motion correction** [[7]](#normcorre) + + * Fast parallelizable OpenCV and FFT-based motion correction of large movies + * Can be run also in online mode (i.e. one frame at a time) + * Corrects for non-rigid artifacts due to raster scanning or non-uniform brain motion + * FFTs can be computed on GPUs (experimental). Requires pycuda and skcuda to be installed. + +* **Source extraction** + + * Separates different sources based on constrained nonnegative matrix Factorization (CNMF) [[1-3]](#caiman) + * Deals with heavily overlapping and neuropil contaminated movies + * Suitable for both 2-photon [[2]](#neuron) and 1-photon [[4]](#cnmfe) calcium imaging data + * Selection of inferred sources using a pre-trained convolutional neural network classifier + * Online processing available [[6]](#onacid) + +* **Denoising, deconvolution and spike extraction** + + * Infers neural activity from fluorescence traces [[2]](#neuron) + * Also works in online mode (i.e. one sample at a time) [[5]](#oasis) + +* **Automatic ROI registration across multiple days** [[1]](#caiman) + +* **Behavioral Analysis** [[8]](#behavior) + + * Unsupervised algorithms based on optical flow and NMF to automatically extract motor kinetics + * Scales to large datasets by exploiting online dictionary learning + * We also developed a tool for acquiring movies at high speed with low cost equipment [[Github repository]](https://github.com/bensondaled/eyeblink). + +* **Variance Stabilization** [[9]](#vst) + * Noise parameters estimation under the Poisson-Gaussian noise model + * Fast algorithm that scales to large datasets + * A basic demo can be found at `CaImAn/demos/notebooks/demo_VST.ipynb` + +# References + +The following references provide the theoretical background and original code for the included methods. + +### Software package detailed description and benchmarking + +If you use this code please cite the corresponding papers where original methods appeared (see References below), as well as: + +[1] Giovannucci A., Friedrich J., Gunn P., Kalfon J., Koay S.A., Taxidis J., Najafi F., Gauthier J.L., Zhou P., Tank D.W., Chklovskii D.B., Pnevmatikakis E.A. (2018). CaImAn: An open source tool for scalable Calcium Imaging data Analysis. bioarXiv preprint. [[paper]](https://doi.org/10.1101/339564) + +### Deconvolution and demixing of calcium imaging data + +[2] Pnevmatikakis, E.A., Soudry, D., Gao, Y., Machado, T., Merel, J., ... & Paninski, L. (2016). Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron 89(2):285-299, [[paper]](http://dx.doi.org/10.1016/j.neuron.2015.11.037), [[Github repository]](https://github.com/epnev/ca_source_extraction). + +[3] Pnevmatikakis, E.A., Gao, Y., Soudry, D., Pfau, D., Lacefield, C., ... & Paninski, L. (2014). A structured matrix factorization framework for large scale calcium imaging data analysis. arXiv preprint arXiv:1409.2903. [[paper]](http://arxiv.org/abs/1409.2903). + +[4] Zhou, P., Resendez, S. L., Stuber, G. D., Kass, R. E., & Paninski, L. (2016). Efficient and accurate extraction of in vivo calcium signals from microendoscopic video data. arXiv preprint arXiv:1605.07266. [[paper]](https://arxiv.org/abs/1605.07266), [[Github repository]](https://github.com/zhoupc/CNMF_E). + +[5] Friedrich J. and Paninski L. Fast active set methods for online spike inference from calcium imaging. NIPS, 29:1984-1992, 2016. [[paper]](https://papers.nips.cc/paper/6505-fast-active-set-methods-for-online-spike-inference-from-calcium-imaging), [[Github repository]](https://github.com/j-friedrich/OASIS). + +### Online Analysis + +[6] Giovannucci, A., Friedrich J., Kaufman M., Churchland A., Chklovskii D., Paninski L., & Pnevmatikakis E.A. (2017). OnACID: Online analysis of calcium imaging data in real data. NIPS 2017, pp. 2378-2388. [[paper]](http://papers.nips.cc/paper/6832-onacid-online-analysis-of-calcium-imaging-data-in-real-time) + +### Motion Correction + +[7] Pnevmatikakis, E.A., and Giovannucci A. (2017). NoRMCorre: An online algorithm for piecewise rigid motion correction of calcium imaging data. Journal of Neuroscience Methods, 291:83-92 [[paper]](https://doi.org/10.1016/j.jneumeth.2017.07.031), [[Github repository]](https://github.com/simonsfoundation/normcorre). + +### Behavioral Analysis + +[8] Giovannucci, A., Pnevmatikakis, E. A., Deverett, B., Pereira, T., Fondriest, J., Brady, M. J., ... & Masip, D. (2017). Automated gesture tracking in head-fixed mice. Journal of Neuroscience Methods, 300:184-195. [[paper]](https://doi.org/10.1016/j.jneumeth.2017.07.014). + +### Variance Stabilization + +[9] Tepper, M., Giovannucci, A., and Pnevmatikakis, E (2018). Anscombe meets Hough: Noise variance stabilization via parametric model estimation. In ICASSP, 2018. [[paper]](https://marianotepper.github.io/papers/anscombe-meets-hough.pdf). [[Github repository]](https://github.com/marianotepper/hough-anscombe) \ No newline at end of file diff --git a/docs/source/CaImAn-features-and-references.md.rst b/docs/source/CaImAn-features-and-references.md.rst new file mode 100644 index 000000000..e226fa568 --- /dev/null +++ b/docs/source/CaImAn-features-and-references.md.rst @@ -0,0 +1,136 @@ +Features +-------- + +CaImAn includes a variety of scalable methods for the analysis of +calcium imaging data: + +- **Handling of very large datasets** + + - Memory mapping + - Parallel processing in patches + - Frame-by-frame online processing `[6] <#onacid>`__ + - OpenCV-based efficient movie playing and resizing + +- **Motion correction** `[7] <#normcorre>`__ + + - Fast parallelizable OpenCV and FFT-based motion correction of + large movies + - Can be run also in online mode (i.e. one frame at a time) + - Corrects for non-rigid artifacts due to raster scanning or + non-uniform brain motion + - FFTs can be computed on GPUs (experimental). Requires pycuda and + skcuda to be installed. + +- **Source extraction** + + - Separates different sources based on constrained nonnegative + matrix Factorization (CNMF) `[1-3] <#caiman>`__ + - Deals with heavily overlapping and neuropil contaminated movies + - Suitable for both 2-photon `[2] <#neuron>`__ and 1-photon + `[4] <#cnmfe>`__ calcium imaging data + - Selection of inferred sources using a pre-trained convolutional + neural network classifier + - Online processing available `[6] <#onacid>`__ + +- **Denoising, deconvolution and spike extraction** + + - Infers neural activity from fluorescence traces `[2] <#neuron>`__ + - Also works in online mode (i.e. one sample at a time) + `[5] <#oasis>`__ + +- **Automatic ROI registration across multiple days** `[1] <#caiman>`__ + +- **Behavioral Analysis** `[8] <#behavior>`__ + + - Unsupervised algorithms based on optical flow and NMF to + automatically extract motor kinetics + - Scales to large datasets by exploiting online dictionary learning + - We also developed a tool for acquiring movies at high speed with + low cost equipment `[Github + repository] `__. + +- **Variance Stabilization** `[9] <#vst>`__ + + - Noise parameters estimation under the Poisson-Gaussian noise model + - Fast algorithm that scales to large datasets + - A basic demo can be found at + ``CaImAn/demos/notebooks/demo_VST.ipynb`` + +References +========== + +The following references provide the theoretical background and original +code for the included methods. + +Software package detailed description and benchmarking +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +If you use this code please cite the corresponding papers where original +methods appeared (see References below), as well as: + +[1] Giovannucci A., Friedrich J., Gunn P., Kalfon J., Koay S.A., Taxidis +J., Najafi F., Gauthier J.L., Zhou P., Tank D.W., Chklovskii D.B., +Pnevmatikakis E.A. (2018). CaImAn: An open source tool for scalable +Calcium Imaging data Analysis. bioarXiv preprint. +`[paper] `__ + +Deconvolution and demixing of calcium imaging data +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +[2] Pnevmatikakis, E.A., Soudry, D., Gao, Y., Machado, T., Merel, J., … +& Paninski, L. (2016). Simultaneous denoising, deconvolution, and +demixing of calcium imaging data. Neuron 89(2):285-299, +`[paper] `__, `[Github +repository] `__. + +[3] Pnevmatikakis, E.A., Gao, Y., Soudry, D., Pfau, D., Lacefield, C., … +& Paninski, L. (2014). A structured matrix factorization framework for +large scale calcium imaging data analysis. arXiv preprint +arXiv:1409.2903. `[paper] `__. + +[4] Zhou, P., Resendez, S. L., Stuber, G. D., Kass, R. E., & Paninski, +L. (2016). Efficient and accurate extraction of in vivo calcium signals +from microendoscopic video data. arXiv preprint arXiv:1605.07266. +`[paper] `__, `[Github +repository] `__. + +[5] Friedrich J. and Paninski L. Fast active set methods for online +spike inference from calcium imaging. NIPS, 29:1984-1992, 2016. +`[paper] `__, +`[Github repository] `__. + +Online Analysis +~~~~~~~~~~~~~~~ + +[6] Giovannucci, A., Friedrich J., Kaufman M., Churchland A., Chklovskii +D., Paninski L., & Pnevmatikakis E.A. (2017). OnACID: Online analysis of +calcium imaging data in real data. NIPS 2017, pp. 2378-2388. +`[paper] `__ + +Motion Correction +~~~~~~~~~~~~~~~~~ + +[7] Pnevmatikakis, E.A., and Giovannucci A. (2017). NoRMCorre: An online +algorithm for piecewise rigid motion correction of calcium imaging data. +Journal of Neuroscience Methods, 291:83-92 +`[paper] `__, `[Github +repository] `__. + +Behavioral Analysis +~~~~~~~~~~~~~~~~~~~ + +[8] Giovannucci, A., Pnevmatikakis, E. A., Deverett, B., Pereira, T., +Fondriest, J., Brady, M. J., … & Masip, D. (2017). Automated gesture +tracking in head-fixed mice. Journal of Neuroscience Methods, +300:184-195. +`[paper] `__. + +Variance Stabilization +~~~~~~~~~~~~~~~~~~~~~~ + +[9] Tepper, M., Giovannucci, A., and Pnevmatikakis, E (2018). Anscombe +meets Hough: Noise variance stabilization via parametric model +estimation. In ICASSP, 2018. +`[paper] `__. +`[Github +repository] `__ diff --git a/docs/source/Getting Started.md b/docs/source/Getting Started.md new file mode 100644 index 000000000..3d8cdd7ee --- /dev/null +++ b/docs/source/Getting Started.md @@ -0,0 +1,72 @@ +## Demos + +* Notebooks: The notebooks provide a simple and friendly way to get into CaImAn and understand its main characteristics. +They are located in the `demos/notebooks`. To launch one of the jupyter notebooks: + + ```bash + source activate CaImAn + jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10 + ``` + and select the notebook from within Jupyter's browser. The argument `--NotebookApp.iopub_data_rate_limit=1.0e10` will prevent any memory issues while plotting on a notebook. + +* demo files are also found in the demos/general subfolder. We suggest trying demo_pipeline.py first as it contains most of the tasks required by calcium imaging. For behavior use demo_behavior.py + +* If you want to directly launch the python files, your python console still must be in the CaImAn directory. + +## Basic Structure + +We recently refactored the code to simplify the parameter setting and usage of the various algorithms. The code now is based revolves around the following objects: + +* `params`: A single object containing a set of dictionaries with the parameters used in all the algorithms. It can be set and changed easily and is passed into all the algorithms. +* `MotionCorrect`: An object for motion correction which can be used for both rigid and piece-wise rigid motion correction. +* `cnmf`: An object for running the CaImAn batch algorithm either in patches or not, suitable for both two-photon (CNMF) and one-photon (CNMF-E) data. +* `online_cnmf`: An object for running the CaImAn online (OnACID) algorithm on two-photon data with or without motion correction. +* `estimates`: A single object that stores the results of the algorithms (CaImAn batch, CaImAn online) in a unified way that also contains plotting methods. + +To see examples of how these methods are used, please consult the demos. + +## Result Interpretation + +As mentioned above, the results of the analysis are stored within the `estimates` objects. The basic entries are the following: + +### Result variables for 2p batch analysis + +The results of CaImAn are saved in an `estimates` object. This is stored inside the cnmf object, i.e. it can be accessed using `cnmf.estimates`. The variables of interest are: +- `estimates.A`: Set of spatial components. Saved as a sparse column format matrix with dimensions (# of pixels X # of components). Each column corresponds to a spatial component. +- `estimates.C`: Set of temporal components. Saved as a numpy array with dimensions (# of components X # of timesteps). Each row corresponds to a background component denoised and deconvolved. +- `estimates.b`: Set of background spatial components (for 2p analysis): Saved as a numpy array with dimensions (# of pixels X # of components). Each column corresponds to a spatial background component. +- `estimates.f`: Set of temporal background components (for 2p analysis). Saved as a numpy array with dimensions (# of background components X # of timesteps). Each row corresponds to a temporal background component. +- `estimates.S`: Deconvolved neural activity (spikes) for each component. Saved as a numpy array with dimensions (# of background components X # of timesteps). Each row corresponds to the deconvolved neural activity for the corresponding component. +- `estimates.YrA`: Set or residual components. Saved as a numpy array with dimensions (# of components X # of timesteps). Each row corresponds to the residual signal after denoising the corresponding component in `estimates.C`. +- `estimates.F_dff`: Set of DF/F normalized temporal components. Saved as a numpy array with dimensions (# of components X # of timesteps). Each row corresponds to the DF/F fluorescence for the corresponding component. + +To view the spatial components, their corresponding vectors need first to be reshaped into 2d images. For example if you want to view the i-th component you can type +``` +import matplotlib.pyplot as plt +plt.figure(); plt.imshow(np.reshape(estimates.A[:,i-1].toarray(), dims, order='F')) +``` +where `dims` is a list or tuple that has the dimensions of the FOV. Similarly if you want to plot the trace for the i-th component you can simply type +``` +plt.figure(); plt.plot(estimates.V[i-1]) +``` +The methods `estimates.plot_contours` and `estimates.view_components` can be used to visualize all the components. + +### Variables for component evaluation + +If you use post-screening to evaluate the quality of the components and remove bad components the results are stored in the lists: +- `idx_components`: List containing the indexes of accepted components. +- `idx_components_bad`: List containing the indexes of rejected components. + +These lists can be used to index the results. For example `estimates.A[:,idx_components]` or `estimates.C[idx_components]` will return the accepted spatial or temporal components, respectively. If you want to view the first accepted component you can type +``` +plt.figure(); plt.imshow(np.reshape(estimates.A[:,idx_components[0]].toarray(), dims, order='F')) +plt.figure(); plt.plot(cnm.estimates.C[idx_components[0]]) +``` + +### Variables for 1p processing (CNMF-E) + +The variables for one photon processing are the same, with an additional variable `estimates.W` for the matrix that is used to compute the background using the ring model, and `estimates.b0` for the baseline value for each pixel. + +### Variables for online processing + +The same `estimates` object is also used for the results of online processing, stored in `onacid.estimates`. diff --git a/docs/source/Getting_Started.rst b/docs/source/Getting_Started.rst new file mode 100644 index 000000000..6554afab1 --- /dev/null +++ b/docs/source/Getting_Started.rst @@ -0,0 +1,137 @@ +Demos +===== + +- Notebooks: The notebooks provide a simple and friendly way to get + into CaImAn and understand its main characteristics. They are located + in the ``demos/notebooks``. To launch one of the jupyter notebooks: + + .. code:: bash + + source activate CaImAn + jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10 + + and select the notebook from within Jupyter’s browser. The argument + ``--NotebookApp.iopub_data_rate_limit=1.0e10`` will prevent any + memory issues while plotting on a notebook. + +- demo files are also found in the demos/general subfolder. We suggest + trying demo_pipeline.py first as it contains most of the tasks + required by calcium imaging. For behavior use demo_behavior.py + +- If you want to directly launch the python files, your python console + still must be in the CaImAn directory. + +Basic Structure +=============== + +We recently refactored the code to simplify the parameter setting and +usage of the various algorithms. The code now is based revolves around +the following objects: + +- ``params``: A single object containing a set of dictionaries with the + parameters used in all the algorithms. It can be set and changed + easily and is passed into all the algorithms. +- ``MotionCorrect``: An object for motion correction which can be used + for both rigid and piece-wise rigid motion correction. +- ``cnmf``: An object for running the CaImAn batch algorithm either in + patches or not, suitable for both two-photon (CNMF) and one-photon + (CNMF-E) data. +- ``online_cnmf``: An object for running the CaImAn online (OnACID) + algorithm on two-photon data with or without motion correction. +- ``estimates``: A single object that stores the results of the + algorithms (CaImAn batch, CaImAn online) in a unified way that also + contains plotting methods. + +To see examples of how these methods are used, please consult the demos. + +Result Interpretation +===================== + +As mentioned above, the results of the analysis are stored within the +``estimates`` objects. The basic entries are the following: + +Result variables for 2p batch analysis +-------------------------------------- + +The results of CaImAn are saved in an ``estimates`` object. This is +stored inside the cnmf object, i.e. it can be accessed using +``cnmf.estimates``. The variables of interest are: + +- ``estimates.A``: Set of spatial components. Saved as a sparse column format matrix with + dimensions (# of pixels X # of components). Each column corresponds to a + spatial component. +- ``estimates.C``: Set of temporal components. Saved as a numpy array with dimensions (# of components X # of timesteps). + Each row corresponds to a background component denoised and deconvolved. +- ``estimates.b``: Set of background spatial components (for 2p + analysis): Saved as a numpy array with dimensions (# of pixels X # of + components). Each column corresponds to a spatial background component. +- ``estimates.f``: Set of temporal background components (for 2p + analysis). Saved as a numpy array with dimensions (# of background + components X # of timesteps). Each row corresponds to a temporal + background component. +- ``estimates.S``: Deconvolved neural activity + (spikes) for each component. Saved as a numpy array with dimensions (# + of background components X # of timesteps). Each row corresponds to the + deconvolved neural activity for the corresponding component. +- ``estimates.YrA``: Set or residual components. Saved as a numpy array + with dimensions (# of components X # of timesteps). Each row corresponds + to the residual signal after denoising the corresponding component in + ``estimates.C``. +- ``estimates.F_dff``: Set of DF/F normalized temporal + components. Saved as a numpy array with dimensions (# of components X # + of timesteps). Each row corresponds to the DF/F fluorescence for the + corresponding component. + +To view the spatial components, their corresponding vectors need first +to be reshaped into 2d images. For example if you want to view the i-th +component you can type + +:: + + import matplotlib.pyplot as plt + plt.figure(); plt.imshow(np.reshape(estimates.A[:,i-1].toarray(), dims, order='F')) + +where ``dims`` is a list or tuple that has the dimensions of the FOV. +Similarly if you want to plot the trace for the i-th component you can +simply type + +:: + + plt.figure(); plt.plot(estimates.V[i-1]) + +The methods ``estimates.plot_contours`` and +``estimates.view_components`` can be used to visualize all the +components. + +Variables for component evaluation +---------------------------------- + +If you use post-screening to evaluate the quality of the components and +remove bad components the results are stored in the lists: - +``idx_components``: List containing the indexes of accepted components. +- ``idx_components_bad``: List containing the indexes of rejected +components. + +These lists can be used to index the results. For example +``estimates.A[:,idx_components]`` or ``estimates.C[idx_components]`` +will return the accepted spatial or temporal components, respectively. +If you want to view the first accepted component you can type + +:: + + plt.figure(); plt.imshow(np.reshape(estimates.A[:,idx_components[0]].toarray(), dims, order='F')) + plt.figure(); plt.plot(cnm.estimates.C[idx_components[0]]) + +Variables for 1p processing (CNMF-E) +------------------------------------ + +The variables for one photon processing are the same, with an additional +variable ``estimates.W`` for the matrix that is used to compute the +background using the ring model, and ``estimates.b0`` for the baseline +value for each pixel. + +Variables for online processing +------------------------------- + +The same ``estimates`` object is also used for the results of online +processing, stored in ``onacid.estimates``. diff --git a/docs/source/INSTALL-windows.rst b/docs/source/INSTALL-windows.rst new file mode 100644 index 000000000..ba75826d5 --- /dev/null +++ b/docs/source/INSTALL-windows.rst @@ -0,0 +1,109 @@ +.. role:: raw-latex(raw) + :format: latex +.. + +.. _install-windows: + +Installation on Windows +======================= + +The Windows installation process differs more widely from installation +on Linux or MacOSX and has different issues you may run into. + +Process +======= + +- Increase the maximum size of your pagefile to 64G or more + (http://www.tomshardware.com/faq/id-2864547/manage-virtual-memory-pagefile-windows.html + ) - The Windows memmap interface is sensitive to the maximum setting + and leaving it at the default can cause errors when processing larger + datasets +- Download and install Anaconda (Python 3.x, not 2.x) + http://docs.continuum.io/anaconda/install. Allow the installer to + modify your PATH variable +- Use Conda to install git (With "conda install git") - use of another + commandline git is acceptable, but may lead to issues depending on + default settings +- Install Microsoft Build Tools for Visual Studio 2017 + https://www.visualstudio.com/downloads/#build-tools-for-visual-studio-2017. + Check the "Build Tools" box, and in the detailed view on the right + check the "C/C++ CLI Tools" component too. The specifics of this + occasionally change as Microsoft changes its products and website; + you may need to go off-script. + +Use the following menu item to launch a anaconda-enabled command prompt: +start>programs>anaconda3>anaconda prompt From that prompt. issue the +following commands (if you wish to use the dev branch, you may switch +branches after the clone): + +.. code:: bash + + git clone https://github.com/flatironinstitute/CaImAn + cd CaImAn + conda env create -f environment.yml -n caiman + conda install -n caiman vs2017_win-64 + +At this point you will want to remove a startup script that visual +studio made for your conda environment that can cause conda to crash +while entering the caiman environment. Use the Windows find-file utility +(under the Start Menu) to look for vs2015\_compiler\_vars.bat and/or +vs2015\_compiler\_vars.bat under your home directory. At least one copy +should show up. Delete the version that has +conda:raw-latex:`\envs`:raw-latex:`\caiman `as part of its location. You +may also want to do a search for keras\_activate.bat under your home +directory, find the one in conda:raw-latex:`\envs`:raw-latex:`\caiman`, +and edit it so KERAS\_BACKEND is set to tensorflow rather than theano. +You may then continue the installation. + +.. code:: bash + + activate caiman + pip install . (OR pip install -e . if you want to develop code) + copy caimanmanager.py .. + conda install numba + cd .. + +Setting up a data directory with caimanmanager +============================================== + +Now that you have stepped out of the caiman source directory, you are +ready to make a data directory with code samples and datasets. You will +not use the source tree directory any more. + +:: + + caimanmanager.py install + +or +``python caimanmanager.py install --inplace`` if you used "pip install +-e ." + +This will place that directory under your home directory in a directory +called caiman\_data. If you have, some of the demos or datafiles may +have changed since your last install, to follow API changes. You can +check to see if they have by doing ``caimanmanager.py check``. If they +have not, you may keep using them. If they have, we recommend moving +your old caiman data directory out of the way (or just remove them if +you have no precious data) and doing a new data install as per above. + +If you prefer to manage this information somewhere else, the +``CAIMAN_DATA`` environment variable can be set to customise it. The +caimanmanager tool and other libraries will respect that. + +Setting up environment variables +================================ + +To make the package work *efficiently* and eliminate "crosstalk" between +different processes, run these commands before launching Python: + +.. code:: bash + + set MKL_NUM_THREADS=1 + set OPENBLAS_NUM_THREADS=1 + set KERAS_BACKEND=tensorflow + +The commands should be run every time you enter the caiman conda +environment. We recommend you save these values inside your environment +so you do not have to repeat this process every time. You can do this by +following the instructions +`here `__. diff --git a/docs/source/Installation.md b/docs/source/Installation.md new file mode 100644 index 000000000..8f5d60c7b --- /dev/null +++ b/docs/source/Installation.md @@ -0,0 +1,145 @@ +## Installation + +Download and install Anaconda or Miniconda (Python 3.x version) + +
+ Installation on Windows + + The Windows installation process differs more widely from installation on Linux or MacOSX and has different issues you may run into. + +### Process + * Increase the maximum size of your pagefile to 64G or more (http://www.tomshardware.com/faq/id-2864547/manage-virtual-memory-pagefile-windows.html ) - The Windows memmap interface is sensitive to the maximum setting and leaving it at the default can cause errors when processing larger datasets + * Download and install Anaconda (Python 3.x) . Allow the installer to modify your PATH variable + * Use Conda to install git (With "conda install git") - use of another commandline git is acceptable, but may lead to issues depending on default settings + * Install Microsoft Build Tools for Visual Studio 2017 . Check the "Build Tools" box, and in the detailed view on the right check the "C/C++ CLI Tools" component too. The specifics of this occasionally change as Microsoft changes its products and website; you may need to go off-script. + +Use the following menu item to launch a anaconda-enabled command prompt: start>programs>anaconda3>anaconda prompt +From that prompt. issue the following commands (if you wish to use the dev branch, you may switch branches after the clone): + + ```bash + git clone https://github.com/flatironinstitute/CaImAn + cd CaImAn + conda env create -f environment.yml -n caiman + conda install -n caiman vs2017_win-64 + ``` + +At this point you will want to remove a startup script that visual studio made for your conda environment that can cause conda to crash while entering the caiman environment. Use the Windows find-file utility (under the Start Menu) to look for vs2015_compiler_vars.bat and/or vs2015_compiler_vars.bat under your home directory. At least one copy should show up. Delete the version that has conda\envs\caiman as part of its location. You may also want to do a search for keras_activate.bat under your home directory, find the one in conda\envs\caiman, and edit it so KERAS_BACKEND is set to tensorflow rather than theano. You may then continue the installation. + + ```bash + activate caiman + pip install . (OR pip install -e . if you want to develop code) + copy caimanmanager.py .. + conda install numba + cd .. + ``` + +### Setting up a data directory with caimanmanager + +Now that you have stepped out of the caiman source directory, you are ready to make a data directory with code samples and datasets. You will not use the source tree directory any more. + + ``` + caimanmanager.py install + ``` + or + ``` + python caimanmanager.py install --inplace + ``` + if you used "pip install -e ." + +This will place that directory under your home directory in a directory called caiman_data. If you have, some of the demos or datafiles may have changed since your last install, to follow API changes. You can check to see if they have by doing `caimanmanager.py check`. If they have not, you may keep using them. If they have, we recommend moving your old caiman data directory out of the way (or just remove them if you have no precious data) and doing a new data install as per above. + +If you prefer to manage this information somewhere else, the `CAIMAN_DATA` environment variable can be set to customise it. The caimanmanager tool and other libraries will respect that. + +### Setting up environment variables + +To make the package work *efficiently* and eliminate "crosstalk" between different processes, run these commands before launching Python: + + ```bash + set MKL_NUM_THREADS=1 + set OPENBLAS_NUM_THREADS=1 + set KERAS_BACKEND=tensorflow + ``` + + The commands should be run every time you enter the caiman conda environment. We recommend you save these values inside your environment so you do not have to repeat this process every time. You can do this by following the instructions [here](https://conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html#saving-environment-variables). + +
+ +
+ Installation on MacOS and Linux + + ```bash + git clone https://github.com/flatironinstitute/CaImAn + cd CaImAn/ + conda env create -f environment.yml -n caiman + source activate caiman + pip install . + ``` + If you want to develop code then replace the last command with + ``` + pip install -e . + ``` + If any of these steps gives you errors do not proceed to the following step without resolving it + + #### known issues + + If you recently upgraded to OSX Mojave you may need to perform the following steps before your first install: + + ``` + xcode-select --install + open /Library/Developer/CommandLineTools/Packages/ + ``` + and install the package file you will find in the folder that pops up + + ### Setting up environment variables + + To make the package work *efficiently* and eliminate "crosstalk" between different processes, run these commands before launching Python (this is for Linux and OSX): + + ```bash + export MKL_NUM_THREADS=1 + export OPENBLAS_NUM_THREADS=1 + export KERAS_BACKEND=tensorflow + ``` + The commands should be run every time before launching python. It is recommended that you save these values inside your environment so you don't have to repeat this process every time. You can do this by following the instructions [here](https://conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html#saving-environment-variables). + + ### Setting up caimanmanager + + Once CaImAn is installed, you may want to get a working directory with code samples and datasets; pip installed a caimanmanager.py command that manages this. If you have not installed Caiman before, you can do + ``` + caimanmanager.py install + ``` + or + ``` + python caimanmanager.py install --inplace + ``` + if you used "pip install -e ." + +This will place that directory under your home directory in a directory called caiman_data. If you have, some of the demos or datafiles may have changed since your last install, to follow API changes. You can check to see if they have by doing `caimanmanager.py check`. If they have not, you may keep using them. If they have, we recommend moving your old caiman data directory out of the way (or just remove them if you have no precious data) and doing a new data install as per above. + +If you prefer to manage this information somewhere else, the `CAIMAN_DATA` environment variable can be set to customise it. The caimanmanager tool and other libraries will respect that. +
+ +## Upgrading + +If you already have CaImAn installed with the pip installer (May 2018 or later), but want to upgrade, please follow the procedure below. If you reinstall CaImAn frequently, you can try skip deleting and recreating your Conda environment. In this case you can do only steps 1, 5, and 7 below to update the code. However, if the environment file has changed since your last update this may lead to you not the latest version. + +From the conda environment you used to install CaImAn: +1. `pip uninstall caiman` +2. remove or rename your ~/caiman_data directory +3. Remove your conda environment: `conda env remove -n NAME_OF_YOUR_ENVIRONMENT` +4. Close and reopen your shell (to clear out the old conda environment) +5. Do a `git pull` from inside your CaImAn folder. +6. Recreate and reenter your conda environment as you did in the [README](https://github.com/flatironinstitute/CaImAn) +7. Do a `pip install .` inside that code checkout +8. Run `caimanmanager.py install` to reinstall the data directory (use `--inplace` if you used the `pip install -e .` during your initial installation). + +- If you used the `pip install -e .` option when installing, then you can try updating by simply doing a `git pull`. Again, this might not lead to the latest version of the code if the environment variables have changed. + +- The same applies if you want to modify some internal function of CaImAn. If you used the `pip install -e .` option then you can directly modify it (that's why it's called developer mode). If you used the `pip install .` option then you will need to `pip uninstall caiman` followed by `pip install .` for your changes to take effect. Depending on the functions you're changing so you might be able to skip this step. + +## Installing additional packages + +CaImAn uses the conda-forge conda channel for installing its required packages. If you want to install new packages into your conda environment for CaImAn, it is important that you not mix conda-forge and the defaults channel; we recommend only using conda-forge. To ensure you're not mixing channels, perform the install (inside your environment) as follows: + ```bash + conda install -c conda-forge --override-channels NEW_PACKAGE_NAME + ``` +You will notice that any packages installed this way will mention, in their listing, that they're from conda-forge, with none of them having a blank origin. If you fail to do this, differences between how packages are built in conda-forge versus the default conda channels may mean that some packages (e.g. OpenCV) stop working despite showing as installed. diff --git a/docs/source/Installation.rst b/docs/source/Installation.rst new file mode 100644 index 000000000..776b8f143 --- /dev/null +++ b/docs/source/Installation.rst @@ -0,0 +1,261 @@ +.. role:: raw-latex(raw) + :format: latex +.. + +Installation +============ + +Download and install Anaconda or Miniconda (Python 3.x version) +http://docs.continuum.io/anaconda/install + +Installation on Windows +------------ + +.. raw:: html + +
+ +The Windows installation process differs more widely from installation +on Linux or MacOSX and has different issues you may run into. + +Process +------- + +- Increase the maximum size of your pagefile to 64G or more + (http://www.tomshardware.com/faq/id-2864547/manage-virtual-memory-pagefile-windows.html + ) - The Windows memmap interface is sensitive to the maximum setting + and leaving it at the default can cause errors when processing larger + datasets + + - Download and install Anaconda (Python 3.x) + http://docs.continuum.io/anaconda/install. Allow the installer to + modify your PATH variable + - Use Conda to install git (With “conda install git”) - use of + another commandline git is acceptable, but may lead to issues + depending on default settings + - Install Microsoft Build Tools for Visual Studio 2017 + https://www.visualstudio.com/downloads/#build-tools-for-visual-studio-2017. + Check the “Build Tools” box, and in the detailed view on the right + check the “C/C++ CLI Tools” component too. The specifics of this + occasionally change as Microsoft changes its products and website; + you may need to go off-script. + +Use the following menu item to launch a anaconda-enabled command prompt: +start>programs>anaconda3>anaconda prompt From that prompt. issue the +following commands (if you wish to use the dev branch, you may switch +branches after the clone): + +.. code:: bash + + git clone https://github.com/flatironinstitute/CaImAn + cd CaImAn + conda env create -f environment.yml -n caiman + conda install -n caiman vs2017_win-64 + +At this point you will want to remove a startup script that visual +studio made for your conda environment that can cause conda to crash +while entering the caiman environment. Use the Windows find-file utility +(under the Start Menu) to look for vs2015_compiler_vars.bat and/or +vs2015_compiler_vars.bat under your home directory. At least one copy +should show up. Delete the version that has +conda:raw-latex:`\envs`:raw-latex:`\caiman `as part of its location. You +may also want to do a search for keras_activate.bat under your home +directory, find the one in conda:raw-latex:`\envs`:raw-latex:`\caiman`, +and edit it so KERAS_BACKEND is set to tensorflow rather than theano. +You may then continue the installation. + +.. code:: bash + + activate caiman + pip install . (OR pip install -e . if you want to develop code) + copy caimanmanager.py .. + conda install numba + cd .. + +Setting up a data directory with caimanmanager +---------------------------------------------- + +Now that you have stepped out of the caiman source directory, you are +ready to make a data directory with code samples and datasets. You will +not use the source tree directory any more. + +:: + + caimanmanager.py install + +or ``python caimanmanager.py install --inplace`` if you used “pip +install -e .” + +This will place that directory under your home directory in a directory +called caiman_data. If you have, some of the demos or datafiles may have +changed since your last install, to follow API changes. You can check to +see if they have by doing ``caimanmanager.py check``. If they have not, +you may keep using them. If they have, we recommend moving your old +caiman data directory out of the way (or just remove them if you have no +precious data) and doing a new data install as per above. + +If you prefer to manage this information somewhere else, the +``CAIMAN_DATA`` environment variable can be set to customise it. The +caimanmanager tool and other libraries will respect that. + +Setting up environment variables +-------------------------------- + +To make the package work *efficiently* and eliminate “crosstalk” between +different processes, run these commands before launching Python: + +.. code:: bash + + set MKL_NUM_THREADS=1 + set OPENBLAS_NUM_THREADS=1 + set KERAS_BACKEND=tensorflow + +The commands should be run every time you enter the caiman conda +environment. We recommend you save these values inside your environment +so you do not have to repeat this process every time. You can do this by +following the instructions +`here `__. + +.. raw:: html + +
+ +Installation on MacOS and Linux +------------- + +.. raw:: html + +
+ +.. code:: bash + + git clone https://github.com/flatironinstitute/CaImAn + cd CaImAn/ + conda env create -f environment.yml -n caiman + source activate caiman + pip install . + +If you want to develop code then replace the last command with +``pip install -e .`` If any of these steps gives you errors do not +proceed to the following step without resolving it + +known issues +~~~~~~~~~~~~ + +If you recently upgraded to OSX Mojave you may need to perform the +following steps before your first install: + +:: + + xcode-select --install + open /Library/Developer/CommandLineTools/Packages/ + +and install the package file you will find in the folder that pops up + +.. _setting-up-environment-variables-1: + +Setting up environment variables +-------------------------------- + +To make the package work *efficiently* and eliminate “crosstalk” between +different processes, run these commands before launching Python (this is +for Linux and OSX): + +.. code:: bash + + export MKL_NUM_THREADS=1 + export OPENBLAS_NUM_THREADS=1 + export KERAS_BACKEND=tensorflow + +The commands should be run every time before launching python. It is +recommended that you save these values inside your environment so you +don’t have to repeat this process every time. You can do this by +following the instructions +`here `__. + +Setting up caimanmanager +------------------------ + +Once CaImAn is installed, you may want to get a working directory with +code samples and datasets; pip installed a caimanmanager.py command that +manages this. If you have not installed Caiman before, you can do +``caimanmanager.py install`` or +``python caimanmanager.py install --inplace`` if you used “pip install +-e .” + +This will place that directory under your home directory in a directory +called caiman_data. If you have, some of the demos or datafiles may have +changed since your last install, to follow API changes. You can check to +see if they have by doing ``caimanmanager.py check``. If they have not, +you may keep using them. If they have, we recommend moving your old +caiman data directory out of the way (or just remove them if you have no +precious data) and doing a new data install as per above. + +If you prefer to manage this information somewhere else, the +``CAIMAN_DATA`` environment variable can be set to customise it. The +caimanmanager tool and other libraries will respect that. + +.. raw:: html + +
+ +Upgrading +========= + +If you already have CaImAn installed with the pip installer (May 2018 or +later), but want to upgrade, please follow the procedure below. If you +reinstall CaImAn frequently, you can try skip deleting and recreating +your Conda environment. In this case you can do only steps 1, 5, and 7 +below to update the code. However, if the environment file has changed +since your last update this may lead to you not the latest version. + +From the conda environment you used to install CaImAn: + +1. ``pip uninstall caiman`` + +2. remove or rename your ~/caiman_data directory + +3. Remove your conda environment: ``conda env remove -n NAME_OF_YOUR_ENVIRONMENT`` + +4. Close and reopen your shell (to clear out the old conda environment) + +5. Do a ``git pull`` from inside your CaImAn folder. + +6. Recreate and reenter your conda environment as you did in the installation instructions + +7. Do a ``pip install .`` inside that code checkout + +8. Run ``caimanmanager.py install`` to reinstall the data directory (use ``--inplace`` if you used the ``pip install -e .`` during your initial installation). + +- If you used the ``pip install -e .`` option when installing, then you + can try updating by simply doing a ``git pull``. Again, this might + not lead to the latest version of the code if the environment + variables have changed. + +- The same applies if you want to modify some internal function of + CaImAn. If you used the ``pip install -e .`` option then you can + directly modify it (that’s why it’s called developer mode). If you + used the ``pip install .`` option then you will need to + ``pip uninstall caiman`` followed by ``pip install .`` for your + changes to take effect. Depending on the functions you’re changing so + you might be able to skip this step. + +Installing additional packages +============================== + +CaImAn uses the conda-forge conda channel for installing its required +packages. If you want to install new packages into your conda +environment for CaImAn, it is important that you not mix conda-forge and +the defaults channel; we recommend only using conda-forge. To ensure +you’re not mixing channels, perform the install (inside your +environment) as follows: + +:: + + conda install -c conda-forge --override-channels NEW_PACKAGE_NAME + +You will notice that any packages installed this way will mention, in +their listing, that they’re from conda-forge, with none of them having a +blank origin. If you fail to do this, differences between how packages +are built in conda-forge versus the default conda channels may mean that +some packages (e.g. OpenCV) stop working despite showing as installed. diff --git a/docs/source/On-file-types-and-sizes.md.rst b/docs/source/On-file-types-and-sizes.md.rst new file mode 100644 index 000000000..af17f5d79 --- /dev/null +++ b/docs/source/On-file-types-and-sizes.md.rst @@ -0,0 +1,57 @@ +CaImAn is designed to perform analysis on datasets saved over a single +or multiple files. However maximum efficiency is achieved when each +dataset is saved as a sequence of files of medium size (1-2GBs). Please +note the following: + +- If you’re using TIFF files make sure that the files are saved in + multipage format. This is particularly important as multipage TIFF + files can be indexed and individual frames can be read without + loading the entire file in memory. On the contrary single page TIFFs + would load the entire file before reading an individual frame. This + can cause significant problems in CaImAn in terms of speed and memory + consumption, as a lot of the parallelization (e.g. during motion + correction) happens by passing the path to a file to multiple + processes each of which will only read and process a small part of + it. **Bear in mind that TIFF files of size 4GB or larger saved + through ImageJ/FIJI are automatically save in single page format and + should be avoided**. If you have such a file you can split into + multiple shorter files through ImageJ/FIJI or through CaImAn using + the following script + +:: + + import numpy as np + import caiman as cm + fname = '' # path to file + m = cm.load(fname) # load the file + T = m.shape(0) # total number of frames for the file + L = 1000 # length of each individual file + fileparts = fname.split(".") + for t in np.arange(0,T,L): + m[t:t+L].save((".").join(fileparts[:-1]) + '_' + str(t//L) + '.' + fileparts[-1]) + +HDF5/H5 files in general do not suffer from this problem. + +- Single frame files should be avoided. The reason is that several + functions, e.g. motion correction, memory mapping, are designed to + work on small sets of frames and in general assume that each file has + more than 1 frames. If your data is saved as a series of single frame + files, you should convert them in a single (or multiple) files. You + can do this by using the following script: + +:: + + import os + import glob + import caiman as cm + fld = '' # path to folder where the data is located + fls = glob.glob(os.path.join(fld,'*.tif')) # change tif to the extension you need + fls.sort() # make sure your files are sorted alphanumerically + m = cm.load_movie_chain(fls) + m.save(os.path.join(fld,'data.tif')) + +If the number of frames is too big, you can split into multiple files as +explained above. Make sure that your files are sorted alphanumerically +before combining them. This can be tricky if your files are initially +saved as ’file1.tif, file2.tif, …, file10.tif`. In this case you can +consult this `page `__. diff --git a/docs/source/Overview.rst b/docs/source/Overview.rst index fe9913e60..ef789cadb 100644 --- a/docs/source/Overview.rst +++ b/docs/source/Overview.rst @@ -1,34 +1,51 @@ Overview ========= -Python translation of Constrained Non-negative Matrix Factorization algorithm for source extraction from calcium imaging data. +.. image:: ../LOGOS/Caiman_logo_FI.png + :width: 200px + :align: right -The code implements a method for simultaneous source extraction and spike inference from large scale calcium imaging movies. The code is suitable for the analysis of somatic imaging data. Implementation for the analysis of dendritic/axonal imaging data will be added in the future. +CaImAn is a Python toolbox for large scale **Ca**\ lcium **Im**\ aging data **An**\ alysis and behavioral analysis. -The algorithm is presented in more detail in +CaImAn implements a set of essential methods required in the analysis pipeline of large scale calcium imaging data. Fast and scalable algorithms are implemented for motion correction, source extraction, spike deconvolution, and component registration across multiple days. It is suitable for both two-photon and one-photon fluorescence microscopy data, and can be run in both batch and online modes. CaImAn also contains some routines for the analysis of behavior from video cameras. A list of features as well as relevant references can be found `here +`_. -Pnevmatikakis, E.A., Soudry, D., Gao, Y., Machado, T., Merel, J., ... & Paninski, L. (2016). Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron, in press, http://dx.doi.org/10.1016/j.neuron.2015.11.037 +Companion paper +-------------- -Pnevmatikakis, E.A., Gao, Y., Soudry, D., Pfau, D., Lacefield, C., ... & Paninski, L. (2014). A structured matrix factorization framework for large scale calcium imaging data analysis. arXiv preprint arXiv:1409.2903. http://arxiv.org/abs/1409.2903 +A paper explaining most of the implementation details and benchmarking can be found `here +`_. -Contributors ------------- +:: + + @article{giovannucci2019caiman, + title={CaImAn: An open source tool for scalable Calcium Imaging data Analysis}, + author={Giovannucci, Andrea and Friedrich, Johannes and Gunn, Pat and Kalfon, Jeremie and Brown, Brandon L and Koay, Sue Ann and Taxidis, Jiannis and Najafi, Farzaneh and Gauthier, Jeffrey L and Zhou, Pengcheng and Khakh, Baljit S and Tank, David W and Chklovskii, Dmitri B and Pnevmatikakis, Eftychios A}, + journal={eLife}, + volume={8}, + pages={e38173}, + year={2019}, + publisher={eLife Sciences Publications Limited} + } -Andrea Giovannucci and -Eftychios Pnevmatikakis -Center for Computational Biology, Simons Foundation, New York, NY +Developers/Contributors +------------ + +CaImAn is being developed at the `Flatiron Institute `_ with numerous contributions from the broader community. The main developers are +* Andrea Giovannucci, University of North Carolina at Chapel Hill, previously at Flatiron Institute +* Eftychios A. Pnevmatikakis, Flatiron Institute +* Johannes Friedrich, Flatiron Institute +* Pat Gunn, Flatiron Institute +A complete list of contributors can be found `here `_. Questions, comments, issues ----------------------------- -Please use the gitter chat room (use the button above) for questions and comments and create an issue for any bugs you might encounter. -Important note ----------------- -The implementation of this package is based on the matlab implementation which can be found [here](https://github.com/epnev/ca_source_extraction). Some of the Matlab features are currently lacking, but will be included in future releases. +Please use our `gitter chat room `_ for questions and comments and create an issue on our `repo page `_ for any bugs you might encounter. License -------- diff --git a/docs/source/Updating-CaImAn.md.rst b/docs/source/Updating-CaImAn.md.rst new file mode 100644 index 000000000..b314efbec --- /dev/null +++ b/docs/source/Updating-CaImAn.md.rst @@ -0,0 +1,28 @@ +If you have you not used the pip installation before (established in May +2018) follow the instructions in the +`README `__ file. + +If you already have CaImAn installed with the pip installer, but want to +upgrade, please follow the procedure below. If you reinstall CaImAn +frequently, you can try skip deleting and recreating your Conda +environment. In this case you can do only steps 1, 5, and 7 below to +update the code. However, if the environment file has changed since your +last update this may lead to you not the latest version. + +From the conda environment you used to install CaImAn: 1. +``pip uninstall caiman`` 2. remove or rename your ~/caiman_data +directory 3. Remove your conda environment: +``conda env remove -n NAME_OF_YOUR_ENVIRONMENT`` 4. Close and reopen +your shell (to clear out the old conda environment) 5. Do a ``git pull`` +from inside your CaImAn folder. 6. Recreate and reenter your conda +environment as you did in the +`README `__ 7. Do a +``pip install .`` inside that code checkout 8. Run +``caimanmanager.py install`` to reinstall the data directory (use +``--inplace`` if you used the ``pip install -e .`` during your initial +installation). + +- If you used the ``pip install -e .`` option, then you can try + updating by simply doing a ``git pull``. Again, this might not lead + to the latest version of the code if the environment variables have + changed. diff --git a/docs/source/conf.py b/docs/source/conf.py index f443ac71b..e1d0bccb5 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -56,7 +56,7 @@ # General information about the project. project = 'CaImAn' -copyright = '2018, Eftychios Pnevmatikakis and Andrea Giovannucci' +copyright = '2019, Flatiron Institute, Simons Foundation, New York, NY' author = 'Eftychios Pnevmatikakis and Andrea Giovannucci' # The version info for the project you're documenting, acts as replacement for @@ -64,9 +64,9 @@ # built documents. # # The short X.Y version. -version = '1.0.0' +version = '1.4' # The full version, including alpha/beta/rc tags. -release = '1.0.0' +release = '1.4' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/source/index.rst b/docs/source/index.rst index 6b56c3ab1..36c2b603f 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -12,6 +12,8 @@ Contents: :maxdepth: 2 Overview + Installation + Getting_Started core_functions diff --git a/environment.yml b/environment.yml index b6c3194f0..338e37249 100644 --- a/environment.yml +++ b/environment.yml @@ -13,6 +13,7 @@ dependencies: - jupyter - keras - matplotlib +- mypy - nose - numpy - opencv diff --git a/mypy.ini b/mypy.ini new file mode 100644 index 000000000..e5a1ca560 --- /dev/null +++ b/mypy.ini @@ -0,0 +1,6 @@ +# See http://mypy.readthedocs.io/en/latest/config_file.html + +[mypy] +check_untyped_defs = True +ignore_missing_imports = True +warn_incomplete_stub = False diff --git a/readthedocs.yml b/readthedocs.yml index 488176203..cae4ae18d 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -1,6 +1,6 @@ # .readthedocs.yml build: - image:latest + image: latest python: version: 3.6 diff --git a/sandbox/demo_caiman_wip.py b/sandbox/demo_caiman_wip.py deleted file mode 100755 index d50390946..000000000 --- a/sandbox/demo_caiman_wip.py +++ /dev/null @@ -1,347 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" -Created on Wed Feb 24 18:39:45 2016 -@author: Andrea Giovannucci -For explanation consult at https://github.com/agiovann/Constrained_NMF/releases/download/v0.4-alpha/Patch_demo.zip -and https://github.com/agiovann/Constrained_NMF -""" -from __future__ import print_function -#%% -from builtins import str -from builtins import range - -try: - if __IPYTHON__: - # this is used for debugging purposes only. allows to reload classes when changed - get_ipython().magic('load_ext autoreload') - get_ipython().magic('autoreload 2') -except NameError: - print('Not launched under iPython') - -import sys -import numpy as np -from time import time -from scipy.sparse import coo_matrix -import psutil -import glob -import os -import scipy -from ipyparallel import Client -import matplotlib as mpl -# mpl.use('Qt5Agg') -import pylab as pl -pl.ion() -#%% -from skimage.external import tifffile -import caiman as cm -from caiman.source_extraction import cnmf as cnmf -from caiman.components_evaluation import evaluate_components -from caiman.utils.visualization import plot_contours -from caiman.base.rois import extract_binary_masks_blob -#%% -fname = '/mnt/ceph/neuro/labeling/k37_20160109_AM_150um_65mW_zoom2p2_00001_1-16/images/k37_20160109_AM_150um_65mW_zoom2p2_00001_00001.tif' -# fname='/run/media/agiovann/ANDREA/k26_v1_176um_target_pursuit_002_013.tif' -fname = './example_movies/demoMovie.tif' - -with tifffile.TiffFile(fname) as tf: - m = cm.movie(tf, fr=30) - mh = cm.movie(tf.asarray(), fr=30) -#(mh-np.percentile(mh,8)).play(backend='opencv',fr=100,gain=10.,magnification=5) -# mh=mh.motion_correct(5,5)[0] -#%% -from sklearn.metrics.pairwise import euclidean_distances -# mh=mh[:,30:60,30:60] -T, h, w = np.shape(mh) -idxA, idxB = np.meshgrid(list(range(w)), list(range(h))) -coordmat = np.vstack((idxA.flatten(), idxB.flatten())) -distanceMatrix = euclidean_distances(coordmat.T) -distanceMatrix = np.logical_and( - distanceMatrix < 2, distanceMatrix > 0) * distanceMatrix -img = np.corrcoef(mh.to_2D(order='C').T) -img = img * distanceMatrix -img = scipy.sparse.coo_matrix(img) -from sklearn.cluster import DBSCAN, KMeans, SpectralClustering -cl = SpectralClustering(n_clusters=3).fit(img) - -#%% -pl.imshow(cl.labels_.reshape((h, w))) -# pl.imshow(mh.local_correlations(eight_neighbours=True,swap_dim=False),alpha=.5) -#%% -pl.imshow(img, interpolation='none') -#%% -cc = scipy.sparse.csgraph.connected_components(img, directed=False) -# distanceMatrix=distanceMatrix/np.max(distanceMatrix) -#%% -m = load(fname, fr=30).resize(1, 1, .2) -#%% -shifts_, xcorrs_, template_ = m.motion_correction_online( - init_frames_template=100, max_shift_h=5, max_shift_w=5) -#%% -shifts, xcorrs, template = m.motion_correction_online( - template=template_, min_count=len(m), max_shift_h=5, max_shift_w=5) -#%% -mov_path = m.apply_shifts_online(shifts, save_base_name='/tmp/test') -#%% -m = load(mov_path, fr=30 * .2) -m1 = cm.concatenate([m.resize(1, 1, 1), mh.resize(1, 1, 1)], axis=1) - -(m1 - np.percentile(m1, 8)).play(backend='opencv', fr=100, gain=3., magnification=5) - -#%% -# mc=m.motion_correct(5,5) -#%% -# m1=mc[0].resize(1,1,.2) -#%% -(m1 - np.percentile(m1, 8)).play(backend='opencv', fr=100, gain=10., magnification=5) -#%% -final_frate = 15 -is_patches = True -is_dendrites = False - - -if is_dendrites == True: - # THIS METHOd CAN GIVE POSSIBLY INCONSISTENT RESULTS ON SOMAS WHEN NOT USED WITH PATCHES - init_method = 'sparse_nmf' - alpha_snmf = None # this controls sparsity -else: - init_method = 'greedy_roi' - alpha_snmf = None - -#%% -# backend='SLURM' -backend = 'local' -if backend == 'SLURM': - n_processes = np.int(os.environ.get('SLURM_NPROCS')) -else: - # roughly number of cores on your machine minus 1 - n_processes = np.maximum(np.int(psutil.cpu_count()), 1) -print(('using ' + str(n_processes) + ' processes')) -#%% start cluster for efficient computation -single_thread = False - -if single_thread: - dview = None -else: - try: - c.close() - except: - print('C was not existing, creating one') - - print("Stopping cluster to avoid unnencessary use of memory....") - sys.stdout.flush() - - if backend == 'SLURM': - try: - stop_server(is_slurm=True) - except: - - print('Nothing to stop') - slurm_script = '/mnt/xfs1/home/agiovann/SOFTWARE/Constrained_NMF/SLURM/slurmStart.sh' - start_server(slurm_script=slurm_script) - - pdir, profile = os.environ['IPPPDIR'], os.environ['IPPPROFILE'] - c = Client(ipython_dir=pdir, profile=profile) - else: - - stop_server() - start_server() - c = Client() - - print(('Using ' + str(len(c)) + ' processes')) - dview = c[:len(c)] -#%% FOR LOADING ALL TIFF FILES IN A FILE AND SAVING THEM ON A SINGLE MEMORY MAPPABLE FILE -fnames = [] -base_folder = '/run/media/agiovann/ANDREA/' # folder containing the demo files -base_folder = '/run/media/agiovann/ANDREA/' -base_folder = './example_movies' -for file in glob.glob(os.path.join(base_folder, '*.tif')): - - if file.endswith("ie.tif"): - fnames.append(os.path.abspath(file)) -fnames.sort() -if len(fnames) == 0: - raise Exception("Could not find any tiff file") - -print(fnames) -fnames = fnames -#%% -# idx_x=slice(12,500,None) -# idx_y=slice(12,500,None) -# idx_xy=(idx_x,idx_y) - -add_to_movie = 300 -downsample_factor = 1. # use .2 or .1 if file is large and you want a quick answer -final_frate = final_frate * downsample_factor -idx_xy = None -base_name = 'Yr' -name_new = save_memmap_each(fnames, dview=dview, base_name=base_name, resize_fact=( - 1, 1, downsample_factor), remove_init=0, idx_xy=idx_xy, add_to_movie=add_to_movie) - -name_new.sort() -print(name_new) -#%% - -fname_new = save_memmap_join( - name_new, base_name='Yr', n_chunks=12, dview=dview) -#%% -# fname_new='Yr_d1_501_d2_398_d3_1_order_F_frames_369_.mmap' -# fname_new=m1 -Yr, dims, T = load_memmap(fname_new) -d1, d2 = dims -images = np.reshape(Yr.T, [T] + list(dims), order='F') -Y = np.reshape(Yr, dims + (T,), order='F') - -#%% -Cn = cnmf.utilities.local_correlations(Y[:, :, :3000]) -pl.imshow(Cn, cmap='gray') - -#%% -if not is_patches: - #%% - - K = 35 # number of neurons expected per patch - gSig = [7, 7] # expected half size of neurons - merge_thresh = 0.8 # merging threshold, max correlation allowed - p = 2 # order of the autoregressive system - cnm = cnmf.CNMF(n_processes, method_init=init_method, k=K, gSig=gSig, merge_thresh=merge_thresh, - p=p, dview=dview, Ain=None) - cnm = cnm.fit(images) - - -#%% -else: - #%% - rf = 15 # half-size of the patches in pixels. rf=25, patches are 50x50 - stride = 4 # amounpl.it of overlap between the patches in pixels - K = 6 # number of neurons expected per patch - gSig = [7, 7] # expected half size of neurons - merge_thresh = 0.8 # merging threshold, max correlation allowed - p = 2 # order of the autoregressive system - memory_fact = 1 # unitless number accounting how much memory should be used. You will need to try different values to see which one would work the default is OK for a 16 GB system - save_results = False - #%% RUN ALGORITHM ON PATCHES - - cnm = cnmf.CNMF(n_processes, k=K, gSig=gSig, merge_thresh=0.8, p=0, dview=c[:], Ain=None, rf=rf, stride=stride, memory_fact=memory_fact, - method_init=init_method, alpha_snmf=alpha_snmf, only_init_patch=True, gnb=1) - cnm = cnm.fit(images) - - A_tot = cnm.A - C_tot = cnm.C - YrA_tot = cnm.YrA - b_tot = cnm.b - f_tot = cnm.f - sn_tot = cnm.sn - - print(('Number of components:' + str(A_tot.shape[-1]))) - - #%% - tB = np.minimum(-2, np.floor(-5. / 30 * final_frate)) - tA = np.maximum(5, np.ceil(25. / 30 * final_frate)) - Npeaks = 10 - traces = C_tot + YrA_tot - # traces_a=traces-scipy.ndimage.percentile_filter(traces,8,size=[1,np.shape(traces)[-1]/5]) - # traces_b=np.diff(traces,axis=1) - fitness_raw, fitness_delta, erfc_raw, erfc_delta, r_values, significant_samples = evaluate_components( - Y, traces, A_tot, C_tot, b_tot, f_tot, remove_baseline=True, N=5, robust_std=False, Athresh=0.1, Npeaks=Npeaks, tB=tB, tA=tA, thresh_C=0.3) - - idx_components_r = np.where(r_values >= .4)[0] - idx_components_raw = np.where(fitness_raw < -20)[0] - idx_components_delta = np.where(fitness_delta < -10)[0] - - idx_components = np.union1d(idx_components_r, idx_components_raw) - idx_components = np.union1d(idx_components, idx_components_delta) - idx_components_bad = np.setdiff1d(list(range(len(traces))), idx_components) - - print(('Keeping ' + str(len(idx_components)) + - ' and discarding ' + str(len(idx_components_bad)))) - #%% - pl.figure() - crd = plot_contours(A_tot.tocsc()[:, idx_components], Cn, thr=0.9) - - #%% - A_tot = A_tot.tocsc()[:, idx_components] - C_tot = C_tot[idx_components] - #%% - save_results = True - if save_results: - np.savez('results_analysis_patch.npz', A_tot=A_tot, C_tot=C_tot, - YrA_tot=YrA_tot, sn_tot=sn_tot, d1=d1, d2=d2, b_tot=b_tot, f=f_tot) - #%% if you have many components this might take long! - pl.figure() - crd = plot_contours(A_tot, Cn, thr=0.9) - #%% - cnm = cnmf.CNMF(n_processes, k=A_tot.shape, gSig=gSig, merge_thresh=merge_thresh, p=p, dview=dview, Ain=A_tot, Cin=C_tot, - f_in=f_tot, rf=None, stride=None) - -# images_1=np.array(images-np.min(images)) - cnm = cnm.fit(images) - -#%% -A, C, b, f, YrA, sn = cnm.A, cnm.C, cnm.b, cnm.f, cnm.YrA, cnm.sn - -#%% -tB = np.minimum(-2, np.floor(-5. / 30 * final_frate)) -tA = np.maximum(5, np.ceil(25. / 30 * final_frate)) -Npeaks = 10 -traces = C + YrA -# traces_a=traces-scipy.ndimage.percentile_filter(traces,8,size=[1,np.shape(traces)[-1]/5]) -# traces_b=np.diff(traces,axis=1) -fitness_raw, fitness_delta, erfc_raw, erfc_delta, r_values, significant_samples = \ - evaluate_components(Y, traces, A, C, b, f, remove_baseline=True, - N=5, robust_std=False, Athresh=0.1, Npeaks=Npeaks, tB=tB, tA=tA, thresh_C=0.3) - - -idx_components_r = np.where(r_values >= .5)[0] -idx_components_raw = np.where(fitness_raw < -40)[0] -idx_components_delta = np.where(fitness_delta < -20)[0] - - -min_radius = gSig[0] - 2 -masks_ws, idx_blobs, idx_non_blobs = extract_binary_masks_blob( - A.tocsc(), min_radius, dims, num_std_threshold=1, - minCircularity=0.5, minInertiaRatio=0.2, minConvexity=.8) - -idx_components = np.union1d(idx_components_r, idx_components_raw) -idx_components = np.union1d(idx_components, idx_components_delta) -idx_blobs = np.intersect1d(idx_components, idx_blobs) -idx_components_bad = np.setdiff1d(list(range(len(traces))), idx_components) - - -print(' ***** ') -print((len(traces))) -print((len(idx_components))) -print((len(idx_blobs))) -#%% -save_results = True -if save_results: - np.savez(os.path.join(os.path.split(fname_new)[0], 'results_analysis.npz'), Cn=Cn, A=A.todense( - ), C=C, b=b, f=f, YrA=YrA, sn=sn, d1=d1, d2=d2, idx_components=idx_components, idx_components_bad=idx_components_bad) - -#%% visualize components -# pl.figure(); -pl.subplot(1, 3, 1) -crd = plot_contours(A.tocsc()[:, idx_components], Cn, thr=0.9) -pl.subplot(1, 3, 2) -crd = plot_contours(A.tocsc()[:, idx_blobs], Cn, thr=0.9) -pl.subplot(1, 3, 3) -crd = plot_contours(A.tocsc()[:, idx_components_bad], Cn, thr=0.9) -#%% -cnmf.utilities.view_patches_bar(Yr, scipy.sparse.coo_matrix(A.tocsc()[ - :, idx_components]), C[idx_components, :], b, f, dims[0], dims[1], YrA=YrA[idx_components, :], img=Cn) -#%% -cnmf.utilities.view_patches_bar(Yr, scipy.sparse.coo_matrix(A.tocsc()[ - :, idx_components_bad]), C[idx_components_bad, :], b, f, dims[0], dims[1], YrA=YrA[idx_components_bad, :], img=Cn) - -#%% STOP CLUSTER -pl.close() -if not single_thread: - c.close() - stop_server() -#%% -stop_server(is_slurm=(backend == 'SLURM')) -log_files = glob.glob('Yr*_LOG_*') - -for log_file in log_files: - os.remove(log_file) diff --git a/use_cases/collaborations/ToliasLab/demo_seeding_nuclear.py b/use_cases/collaborations/ToliasLab/demo_seeding_nuclear.py index 405d3c3bc..37c9f6289 100644 --- a/use_cases/collaborations/ToliasLab/demo_seeding_nuclear.py +++ b/use_cases/collaborations/ToliasLab/demo_seeding_nuclear.py @@ -2,15 +2,16 @@ # -*- coding: utf-8 -*- """ -Demo for detecting ROIs in a structural channel and then seeding CNMF with them. -Detection happens through simple adaptive thresholding of the mean image and -could potentially be improved. Then the structural channel is processed. Both -offline and online approaches are included. +Demo for detecting ROIs in a structural channel and then seeding CNMF with +them. Detection happens through simple adaptive thresholding of the mean image +and could potentially be improved. Then the structural channel is processed. +Both offline and online approaches are included. -The offline approach will only use the seeded masks, whereas online will also +The offline approach will only use the seeded masks, whereas online will also search for new components during the analysis. -The demo assumes that both channels are motion corrected prior to the analysis. +The demo assumes that both channels are motion corrected prior to the analysis +although this is not necessary for the case of caiman online. @author: epnevmatikakis """ @@ -18,7 +19,6 @@ try: if __IPYTHON__: print('Debugging!') - # this is used for debugging purposes only. allows to reload classes when changed get_ipython().magic('load_ext autoreload') get_ipython().magic('autoreload 2') except NameError: @@ -26,193 +26,176 @@ import numpy as np import glob -import pylab as pl +import matplotlib.pyplot as plt import caiman as cm -from caiman.components_evaluation import evaluate_components -from caiman.source_extraction.cnmf import cnmf as cnmf -from caiman.source_extraction.cnmf.online_cnmf import seeded_initialization +# from caiman.motion_correction import MotionCorrect +from caiman.source_extraction import cnmf as cnmf import os -from copy import deepcopy -from caiman.summary_images import max_correlation_image -#%% construct the seeding matrix using the structural channel (note that some components are missed - thresholding can be improved) +# %% use this cell for splitting channels and performing motion correction +# # %% extract individual channels and motion correct them +# +# fname = '/nuclear/gmc_960_30mw_00001.tif' +# fname_green = fname.split('.')[0] + '_green_raw.tif' +# fname_red = fname.split('.')[0] + '_red_raw.tif' +# m = cm.load(fname) +# +# m[::2].save(fname_green) +# m[1::2].save(fname_red) +# +# +# # %% motion correction +# +# # dataset dependent parameters +# fr = 30 # imaging rate in frames per second +# dxy = (1., 1.) # 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 pathes (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': fname_red, +# 'fr': fr, +# '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) +# +# #%% +# +# c, dview, n_processes = cm.cluster.setup_cluster( +# backend='local', n_processes=None, single_thread=False) +# +# # %%% MOTION CORRECTION +# # first we create a motion correction object with the specified parameters +# mc = MotionCorrect(fname_red, dview=dview, **opts.get_group('motion')) +# # note that the file is not loaded in memory +# +# # %% Run (piecewise-rigid motion) correction using NoRMCorre +# mc.motion_correct(save_movie=True) +# R = cm.load(mc.mmap_file, in_memory=True) +# R.save(fname.split('.')[0] + '_red.tif') +# G = mc.apply_shifts_movie(fname_green) + +# %% construct the seeding matrix using the structural channel (note that some +# components are missed - thresholding can be improved) + +fname_red = ('/Users/epnevmatikakis/Documents/Ca_datasets' + + '/Tolias/nuclear/gmc_960_30mw_00001_red.tif') - -filename = 'example_movies/gmc_960_30mw_00001_red.tif' Ain, mR = cm.base.rois.extract_binary_masks_from_structural_channel( - cm.load(filename), expand_method='dilation', selem=np.ones((1, 1))) -pl.figure() + cm.load(fname_red), expand_method='dilation', selem=np.ones((1, 1))) +plt.figure() crd = cm.utils.visualization.plot_contours( Ain.astype('float32'), mR, thr=0.99, display_numbers=False) -pl.title('Contour plots of detected ROIs in the structural channel') +plt.title('Contour plots of detected ROIs in the structural channel') -#%% choose whether to use online algorithm (OnACID) or offline (CNMF) +# %% choose whether to use online algorithm (OnACID) or offline (CNMF) use_online = True -#%% some common parameters +# specify some common parameters + +fname_green = ('/Users/epnevmatikakis/Documents/Ca_datasets/Tolias' + + '/nuclear/gmc_960_30mw_00001_green.tif') + +# %% some common parameters K = 5 # number of neurons expected per patch (nuisance parameter in this case) gSig = [7, 7] # expected half size of neurons +fr = 30 +decay_time = 0.5 merge_thresh = 0.8 # merging threshold, max correlation allowed p = 1 # order of the autoregressive system -#%% -if use_online: - #%% prepare parameters - fnames = 'example_movies/gmc_980_30mw_00001_green.tif' - rval_thr = .95 - thresh_fitness_delta = -30 - thresh_fitness_raw = -30 - initbatch = 100 # use the first initbatch frames to initialize OnACID - # length of dataset (currently needed to allocate matrices) - T1 = 2000 - expected_comps = 500 # maximum number of components - Y = cm.load(fnames, subindices=slice( - 0, initbatch, None)).astype(np.float32) - Yr = Y.transpose(1, 2, 0).reshape((np.prod(Y.shape[1:]), -1), order='F') - - #%% run seeded initialization - cnm_init = seeded_initialization(Y[:initbatch].transpose(1, 2, 0), Ain=Ain, init_batch=initbatch, gnb=1, - gSig=gSig, merge_thresh=0.8, - p=1, minibatch_shape=100, minibatch_suff_stat=5, - update_num_comps=True, rval_thr=rval_thr, - thresh_fitness_delta=thresh_fitness_delta, - thresh_fitness_raw=thresh_fitness_raw, - batch_update_suff_stat=True, max_comp_update_shape=5, simultaneously=True) - - Cn_init = Y.local_correlations(swap_dim=False) - pl.figure() - crd = cm.utils.visualization.plot_contours( - cnm_init.A.tocsc(), Cn_init, thr=0.9) - pl.title('Contour plots of detected ROIs in the structural channel') - - # contour plot after seeded initialization. Note how the contours are not clean since there is no activity - # for most of the ROIs during the first initbatch frames - #%% run OnACID - cnm = deepcopy(cnm_init) - # prepare the object to run OnACID - cnm._prepare_object(np.asarray(Yr), T1, expected_comps) - cnm.max_comp_update_shape = np.inf - cnm.update_num_comps = True - t = cnm.initbatch - - Y_ = cm.load(fnames, subindices=slice(t, T1, None)).astype(np.float32) - - Cn = max_correlation_image(Y_, swap_dim=False) - - for frame_count, frame in enumerate(Y_): - if frame_count % 100 == 99: - print([frame_count, cnm.Ab.shape]) - -# no motion correction here -# templ = cnm.Ab.dot(cnm.C_on[:cnm.M, t - 1]).reshape(cnm.dims, order='F') -# frame_cor, shift = motion_correct_iteration_fast(frame, templ, max_shift, max_shift) -# shifts.append(shift) - cnm.fit_next(t, frame.copy().reshape(-1, order='F')) - t += 1 - - C = cnm.C_on[cnm.gnb:cnm.M] - A = cnm.Ab[:, cnm.gnb:cnm.M] - print(('Number of components:' + str(A.shape[-1]))) - - #%% plot some results - pl.figure() - crd = cm.utils.visualization.plot_contours(A, Cn, thr=0.9) - pl.title('Contour plots of components in the functional channel') - - #%% view components. Last components are the components added by OnACID - dims = Y.shape[1:] - cm.utils.visualization.view_patches_bar(Yr, A, C, cnm.b, cnm.C_on[:cnm.gnb], - dims[0], dims[1], YrA=cnm.noisyC[cnm.gnb:cnm.M] - C, img=Cn) -#%% -else: # run offline CNMF algorithm - #%% start cluster +gnb = 1 # order of background +min_SNR = 2 # trace SNR threshold +rval_thr = .95 + +# %% create and fit online object + +if use_online: # run OnAcid algorithm (CaImAn online) + + show_movie = not True + init_batch = 200 # use the first initbatch frames to initialize OnACID + + params_dict = {'fnames': fname_green, + 'fr': fr, + 'decay_time': decay_time, + 'gSig': gSig, + 'p': p, + 'min_SNR': min_SNR, + 'rval_thr': rval_thr, + 'nb': gnb, + 'motion_correct': False, + 'init_batch': init_batch, + 'init_method': 'seeded', + 'normalize': True, + 'K': K, + 'dist_shape_update': True, + 'show_movie': show_movie} + opts = cnmf.params.CNMFParams(params_dict=params_dict) + + cnm_on = cnmf.online_cnmf.OnACID(params=opts) + cnm_on.estimates.A = Ain + cnm_on.fit_online() + + # %% plot some results + + cnm_on.estimates.plot_contours() + cnm_on.estimates.view_components() + # view components. Last components are the components added by OnACID + +else: # run offline CNMF algorithm (CaImAn Batch) + # %% start cluster c, dview, n_processes = cm.cluster.setup_cluster( backend='local', n_processes=None, single_thread=False) - #%% FOR LOADING ALL TIFF FILES IN A FILE AND SAVING THEM ON A SINGLE MEMORY MAPPABLE FILE + # %% set parameters, cnmf object, and fit - # can actually be a lost of movie to concatenate - fnames = ['example_movies/gmc_980_30mw_00001_green.tif'] - add_to_movie = 0 # the movie must be positive!!! - downsample_factor = .5 # use .2 or .1 if file is large and you want a quick answer - base_name = 'Yr' - name_new = cm.save_memmap_each(fnames, dview=dview, base_name=base_name, resize_fact=( - 1, 1, downsample_factor), add_to_movie=add_to_movie) - name_new.sort() - fname_new = cm.save_memmap_join(name_new, base_name='Yr', dview=dview) + params_dict = {'fnames': fname_green, + 'fr': fr, + 'decay_time': decay_time, + 'gSig': gSig, + 'p': p, + 'min_SNR': min_SNR, + 'rval_thr': rval_thr, + 'nb': gnb, + 'only_init': False, + 'rf': None} - #%% LOAD MEMORY MAPPABLE FILE + opts = cnmf.params.CNMFParams(params_dict=params_dict) + cnm_b = cnmf.CNMF(n_processes, params=opts, dview=dview, Ain=Ain) + cnm_b.fit_file(motion_correct=False) - Yr, dims, T = cm.load_memmap(fname_new) - d1, d2 = dims + # %% load memory mapped file and evaluate components + + Yr, dims, T = cm.load_memmap(cnm_b.mmap_file) images = np.reshape(Yr.T, [T] + list(dims), order='F') - Y = np.reshape(Yr, dims + (T,), order='F') - - #%% play movie, press q to quit - - play_movie = False - if play_movie: - cm.movie(images).play(fr=50, magnification=3, gain=2.) - - #%% movie cannot be negative! - - if np.min(images) < 0: - raise Exception('Movie too negative, add_to_movie should be larger') - - #%% correlation image. From here infer neuron size and density - - Cn = cm.movie(images)[:3000].local_correlations(swap_dim=False) - - #%% run seeded CNMF - - cnm = cnmf.CNMF(n_processes, method_init='greedy_roi', k=Ain.shape[1], gSig=gSig, merge_thresh=merge_thresh, - p=p, dview=dview, Ain=Ain, method_deconvolution='oasis', rolling_sum=False, rf=None) - cnm = cnm.fit(images) - A, C, b, f, YrA, sn = cnm.A, cnm.C, cnm.b, cnm.f, cnm.YrA, cnm.sn - - #%% plot contours of components - - pl.figure() - crd = cm.utils.visualization.plot_contours(cnm.A, Cn, thr=0.9) - pl.title('Contour plots against correlation image') - - #%% evaluate the quality of the components - # a lot of components will be removed because presumably they are not active - # during these 2000 frames of the experiment - - final_frate = 15 # approximate frame rate of data - Npeaks = 10 - traces = C + YrA - fitness_raw, fitness_delta, erfc_raw, erfc_delta, r_values, significant_samples = \ - evaluate_components(Y, traces, A, C, b, f, final_frate, remove_baseline=True, - N=5, robust_std=False, Npeaks=Npeaks, thresh_C=0.3) - - # filter based on spatial consistency - idx_components_r = np.where(r_values >= .85)[0] - # filter based on transient size - idx_components_raw = np.where(fitness_raw < -50)[0] - # filter based on transient derivative size - idx_components_delta = np.where(fitness_delta < -50)[0] - idx_components = np.union1d(idx_components_r, idx_components_raw) - idx_components = np.union1d(idx_components, idx_components_delta) - idx_components_bad = np.setdiff1d(list(range(len(traces))), idx_components) - print((len(traces))) - print((len(idx_components))) - - #%% visualize components - # pl.figure(); - pl.subplot(1, 2, 1) - crd = cm.utils.visualization.plot_contours( - A.tocsc()[:, idx_components], Cn, thr=0.9) - pl.title('selected components') - pl.subplot(1, 2, 2) - pl.title('discaded components') - crd = cm.utils.visualization.plot_contours( - A.tocsc()[:, idx_components_bad], Cn, thr=0.9) - #%% visualize selected components - cm.utils.visualization.view_patches_bar(Yr, A.tocsc()[:, idx_components], C[ - idx_components, :], b, f, dims[0], dims[1], YrA=YrA[idx_components, :], img=Cn) - - #%% STOP CLUSTER and clean up log files + Cn = cm.local_correlations(images, swap_dim=False) + cnm_b.estimates.plot_contours(img=Cn) + + # %% evalute components and do some plotting + cnm_b.estimates.evaluate_components(images, cnm_b.params, dview=dview) + cnm_b.estimates.plot_contours(img=Cn, idx=cnm_b.estimates.idx_components) + cnm_b.estimates.view_components(images, img=Cn, + idx=cnm_b.estimates.idx_components) + + # %% STOP CLUSTER and clean up log files cm.stop_server() log_files = glob.glob('Yr*_LOG_*')