diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml index 52325b3dd..c33caedaf 100644 --- a/.github/workflows/run_tests.yml +++ b/.github/workflows/run_tests.yml @@ -30,6 +30,8 @@ jobs: activate-environment: caiman conda-solver: libmamba miniforge-version: latest + channels: conda-forge + conda-remove-defaults: "true" - name: Install OS Dependencies shell: bash -l {0} diff --git a/README.md b/README.md index b7141dc8e..afbd9e0e0 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ There are two primary ways to install Caiman. The easiest route is to install the miniforge distribution of Anaconda, and use that to install the rest using prebuilt packages. Most users should take this path. ## Route B -The alternative route is to make sure you have a working compiler, create a python virtualenv, grab the caiman sources, and use pip to populate the virtualenv and build Caiman. This route is not as tested and is not presently documented; it is a standard pip-based install. +The alternative route is to make sure you have a working compiler, create a python virtualenv, grab the caiman sources, and use pip to populate the virtualenv and build Caiman. This route is not as tested and is not presently documented; it is a standard pip-based install (although it will invoke your C++ compiler to build some components). # Quick start (Route A) Follow these three steps to get started quickly, from installation to working through a demo notebook. If you do not already have conda installed, [you can find it here](https://github.com/conda-forge/miniforge). The miniforge distribution of conda is preferred; it will require fewer steps and likely encounter fewer issues. If you are using a different distro of conda, you will likely need to add `-c conda-forge` to the commands you use to make your environment. @@ -62,6 +62,11 @@ Jupyter will open. Navigate to demos/notebooks/ and click on `demo_pipeline.ipyn > `` in the first line is your home directory, its location depdnding on your OS/computer. On Linux/Mac it is `~` while on Windows it will be something like `C:\Users\your_user_name\` +# Quick Start (Route B) +This differs from the quick start above in two ways: +* For the first step only, go to [this doc](https://github.com/flatironinstitute/CaImAn/blob/main/docs/source/Installation.rst) and run through the parts of section 1B relevant to your operating system. After that, steps 2 and onward are the same +* You will probably want to manually set some environment variables before any use of caiman; see [here](https://github.com/conda-forge/caiman-feedstock/blob/main/recipe/activate.sh) for a Linux/OSX example, or [here](https://github.com/conda-forge/caiman-feedstock/blob/main/recipe/activate.bat) for a Windows example. Either make a note of this or modify your dotfiles/configuration to do it for you. + ## For installation help Caiman should install easily on Linux, Mac, and Windows. If you run into problems, we have a dedicated [installation page](./docs/source/Installation.rst). If you don't find what you need there, [create an issue](https://github.com/flatironinstitute/Caiman/issues) on GitHub. diff --git a/VERSION b/VERSION index 3d0e62313..5a7dde332 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.11.4 +1.12.1+dj1 diff --git a/caiman/base/movies.py b/caiman/base/movies.py index f0a9a5796..ccf71e73b 100644 --- a/caiman/base/movies.py +++ b/caiman/base/movies.py @@ -172,90 +172,6 @@ def motion_correct(self, return self, shifts, xcorrs, template - def motion_correct_3d(self, - max_shift_z=5, - max_shift_w=5, - max_shift_h=5, - num_frames_template=None, - template=None, - method='opencv', - remove_blanks=False, - interpolation='cubic'): - """ - Extract shifts and motion corrected movie automatically, - - for more control consider the functions extract_shifts and apply_shifts - Disclaimer, it might change the object itself. - - Args: - max_shift,z,max_shift_w,max_shift_h: maximum pixel shifts allowed when - correcting in the axial, width, and height directions - - template: if a good template for frame by frame correlation exists - it can be passed. If None it is automatically computed - - method: depends on what is installed 'opencv' or 'skimage'. 'skimage' - is an order of magnitude slower - - num_frames_template: if only a subset of the movies needs to be loaded - for efficiency/speed reasons - - - Returns: - self: motion corrected movie, it might change the object itself - - shifts : tuple, contains x, y, and z shifts and correlation with template - - xcorrs: cross correlation of the movies with the template - - template: the computed template - """ - - if template is None: # if template is not provided it is created - if num_frames_template is None: - num_frames_template = 10e7 / (self.shape[1] * self.shape[2]) - - frames_to_skip = int(np.maximum(1, self.shape[0] / num_frames_template)) - - # sometimes it is convenient to only consider a subset of the - # movie when computing the median - submov = self[::frames_to_skip, :].copy() - templ = submov.bin_median_3d() # create template with portion of movie - shifts, xcorrs = submov.extract_shifts_3d( - max_shift_z=max_shift_z, # NOTE: extract_shifts_3d has not been implemented yet - use skimage - max_shift_w=max_shift_w, - max_shift_h=max_shift_h, - template=templ, - method=method) - submov.apply_shifts_3d( # NOTE: apply_shifts_3d has not been implemented yet - shifts, interpolation=interpolation, method=method) - template = submov.bin_median_3d() - del submov - m = self.copy() - shifts, xcorrs = m.extract_shifts_3d(max_shift_z=max_shift_z, - max_shift_w=max_shift_w, - max_shift_h=max_shift_h, - template=template, - method=method) - m = m.apply_shifts_3d(shifts, interpolation=interpolation, method=method) - template = (m.bin_median_3d()) - del m - else: - template = template - np.percentile(template, 8) - - # now use the good template to correct - shifts, xcorrs = self.extract_shifts_3d(max_shift_z=max_shift_z, - max_shift_w=max_shift_w, - max_shift_h=max_shift_h, - template=template, - method=method) - self = self.apply_shifts_3d(shifts, interpolation=interpolation, method=method) - - if remove_blanks: - raise Exception("motion_correct_3d(): The remove_blanks parameter was never functional and should not be used") - - return self, shifts, xcorrs, template - def bin_median(self, window: int = 10) -> np.ndarray: """ compute median of 3D array in along axis o by binning values @@ -465,7 +381,6 @@ def apply_shifts(self, shifts, interpolation: str = 'linear', method: str = 'ope min_, max_) elif method == 'skimage': - tform = skimage.transform.AffineTransform(translation=(-sh_y_n, -sh_x_n)) self[i] = skimage.transform.warp(frame, tform, preserve_range=True, order=interpolation) @@ -477,54 +392,6 @@ def apply_shifts(self, shifts, interpolation: str = 'linear', method: str = 'ope return self - def debleach(self): - """ Debleach by fiting a model to the median intensity. - """ - #todo: todocument - if not isinstance(self[0, 0, 0], np.float32): - warnings.warn('Casting the array to float32') - self = np.asanyarray(self, dtype=np.float32) - - t, _, _ = self.shape - x = np.arange(t) - y = np.median(self.reshape(t, -1), axis=1) - - def expf(x, a, b, c): - return a * np.exp(-b * x) + c - - def linf(x, a, b): - return a * x + b - - try: - 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: - p0 = (float(y[-1] - y[0]) / float(x[-1] - x[0]), y[0]) - popt, _ = scipy.optimize.curve_fit(linf, x, y, p0=p0) - y_fit = linf(x, *popt) - - norm = y_fit - np.median(y[:]) - for frame in range(t): - self[frame, :, :] = self[frame, :, :] - norm[frame] - - return self - - def return_cropped(self, crop_top=0, crop_bottom=0, crop_left=0, crop_right=0, crop_begin=0, crop_end=0) -> np.ndarray: - """ - Return a cropped version of the movie - The returned version is independent of the original, which is less memory-efficient but also less likely to be surprising. - - Args: - crop_top/crop_bottom/crop_left,crop_right: how much to trim from each side - - crop_begin/crop_end: (undocumented) - """ - t, h, w = self.shape - ret = np.zeros(( t - crop_end - crop_begin, h - crop_bottom - crop_top, w - crop_right - crop_left)) - ret[:,:,:] = self[crop_begin:t - crop_end, crop_top:h - crop_bottom, crop_left:w - crop_right] - return ret - def removeBL(self, windowSize:int=100, quantilMin:int=8, in_place:bool=False, returnBL:bool=False): """ Remove baseline from movie using percentiles over a window @@ -655,30 +522,6 @@ def computeDFF(self, secsWindow: int = 5, quantilMin: int = 8, method: str = 'on meta_data=self.meta_data, file_name=self.file_name) - def NonnegativeMatrixFactorization(self, - n_components: int = 30, - init: str = 'nndsvd', - beta: int = 1, - tol=5e-7, - sparseness: str = 'components', - **kwargs) -> tuple[np.ndarray, np.ndarray]: - """ - See documentation for scikit-learn NMF - """ - if np.min(self) < 0: - raise ValueError("All values must be positive") - - T, h, w = self.shape - Y = np.reshape(self, (T, h * w)) - Y = Y - np.percentile(Y, 1) - Y = np.clip(Y, 0, np.Inf) - estimator = sklearn.decomposition.NMF(n_components=n_components, init=init, tol=tol, **kwargs) - time_components = estimator.fit_transform(Y) - components_ = estimator.components_ - space_components = np.reshape(components_, (n_components, h, w)) - - return space_components, time_components - def IPCA(self, components: int = 50, batch: int = 1000) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ Iterative Principal Component analysis, see sklearn.decomposition.incremental_pca @@ -855,46 +698,6 @@ def local_correlations(self, return Cn - def partition_FOV_KMeans(self, - tradeoff_weight: float = .5, - fx: float = .25, - fy: float = .25, - n_clusters: int = 4, - max_iter: int = 500) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - Partition the FOV in clusters that are grouping pixels close in space and in mutual correlation - - Args: - tradeoff_weight:between 0 and 1 will weight the contributions of distance and correlation in the overall metric - - fx,fy: downsampling factor to apply to the movie - - n_clusters,max_iter: KMeans algorithm parameters - - Returns: - fovs:array 2D encoding the partitions of the FOV - - mcoef: matrix of pairwise correlation coefficients - - distanceMatrix: matrix of picel distances - """ - _, h1, w1 = self.shape - self.resize(fx, fy) - T, h, w = self.shape - Y = np.reshape(self, (T, h * w)) - mcoef = np.corrcoef(Y.T) - - idxA, idxB = np.meshgrid(list(range(w)), list(range(h))) - coordmat = np.vstack((idxA.flatten(), idxB.flatten())) - distanceMatrix = sklearn.metrics.pairwise.euclidean_distances(coordmat.T) - distanceMatrix = distanceMatrix / np.max(distanceMatrix) - estim = sklearn.cluster.KMeans(n_clusters=n_clusters, max_iter=max_iter) - kk = estim.fit(tradeoff_weight * mcoef - (1 - tradeoff_weight) * distanceMatrix) - labs = kk.labels_ - fovs = np.reshape(labs, (h, w)) - fovs = cv2.resize(np.uint8(fovs), (w1, h1), 1. / fx, 1. / fy, interpolation=cv2.INTER_NEAREST) - return np.uint8(fovs), mcoef, distanceMatrix - def extract_traces_from_masks(self, masks: np.ndarray) -> caiman.base.traces.trace: """ Args: @@ -978,19 +781,6 @@ def resize(self, fx=1, fy=1, fz=1, interpolation=cv2.INTER_AREA): return self - def guided_filter_blur_2D(self, guide_filter, radius: int = 5, eps=0): - """ - performs guided filtering on each frame. See opencv documentation of cv2.ximgproc.guidedFilter - """ - logger = logging.getLogger("caiman") - - for idx, fr in enumerate(self): - if idx % 1000 == 0: - logger.debug(f"At index: {idx}") - self[idx] = cv2.ximgproc.guidedFilter(guide_filter, fr, radius=radius, eps=eps) - - return self - def bilateral_blur_2D(self, diameter: int = 5, sigmaColor: int = 10000, sigmaSpace=0): """ performs bilateral filtering on each frame. See opencv documentation of cv2.bilateralFilter @@ -1071,33 +861,6 @@ def to_2D(self, order='F') -> np.ndarray: T = self.shape[0] return np.reshape(self, (T, -1), order=order) - def zproject(self, method: str = 'mean', cmap=matplotlib.cm.gray, aspect='auto', **kwargs) -> np.ndarray: - """ - Compute and plot projection across time: - - Args: - method: String - 'mean','median','std' - - **kwargs: dict - arguments to imagesc - - Raises: - Exception 'Method not implemented' - """ - # TODO: make the imshow optional - # TODO: todocument - if method == 'mean': - zp = np.mean(self, axis=0) - elif method == 'median': - zp = np.median(self, axis=0) - elif method == 'std': - zp = np.std(self, axis=0) - else: - raise Exception('Method not implemented') - plt.imshow(zp, cmap=cmap, aspect=aspect, **kwargs) - return zp - def play(self, gain: float = 1, fr=None, @@ -2340,3 +2103,4 @@ def view(button): def rgb2gray(rgb): # Standard mathematical conversion return np.dot(rgb[..., :3], [0.299, 0.587, 0.114]) + diff --git a/caiman/caimanmanager.py b/caiman/caimanmanager.py index 968bbfca6..d0e00969a 100755 --- a/caiman/caimanmanager.py +++ b/caiman/caimanmanager.py @@ -45,7 +45,8 @@ # standard_movies: These are needed by the demo standard_movies = [ os.path.join('example_movies', 'data_endoscope.tif'), - os.path.join('example_movies', 'demoMovie.tif') + os.path.join('example_movies', 'demoMovie.tif'), + os.path.join('example_movies', 'avg_mask_fixed.png') ] ############### diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py index 9b3a4e7fb..f5c030dbf 100644 --- a/caiman/motion_correction.py +++ b/caiman/motion_correction.py @@ -472,6 +472,7 @@ def apply_shifts_movie(self, fname, rigid_shifts: Optional[bool] = None, save_me if save_memmap: dims = m_reg.shape fname_tot = caiman.paths.memmap_frames_filename(save_base_name, dims[1:], dims[0], order) + fname_tot = caiman.paths.fn_relocated(fname_tot) big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=caiman.mmapping.prepare_shape((np.prod(dims[1:]), dims[0])), order=order) big_mov[:] = np.reshape(m_reg.transpose(1, 2, 0), (np.prod(dims[1:]), dims[0]), order='F') diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index c8f8d6bdc..8c9c5ec58 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -290,7 +290,32 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p self.estimates = Estimates(A=Ain, C=Cin, b=b_in, f=f_in, dims=self.params.data['dims']) - def fit_file(self, motion_correct=False, indices=None, include_eval=False): + def __str__(self): + ret = f"Caiman CNMF Object. subfields:{list(self.__dict__.keys()) }" + if hasattr(self.estimates, 'A') and self.estimates.A is not None: + ret += f" A.shape={self.estimates.A.shape}" + if hasattr(self.estimates, 'b') and self.estimates.b is not None: + ret += f" b.shape={self.estimates.b.shape}" + if hasattr(self.estimates, 'C') and self.estimates.C is not None: + ret += f" C.shape={self.estimates.C.shape}" + return ret + + def __repr__(self): + ret = f"Caiman CNMF Object" + if hasattr(self.estimates, 'A') and self.estimates.A is not None: + ret += f" A.shape={self.estimates.A.shape}" + if hasattr(self.estimates, 'b') and self.estimates.b is not None and len(self.estimates.b.shape) > 1: + ret += f" bg components={self.estimates.b.shape[1]}" + if hasattr(self.estimates, 'C') and self.estimates.C is not None: + ret += f" C.shape={self.estimates.C.shape}" + ret += " Use str() for more details" + return ret + + def __getitem__(self, idx): + return getattr(self, idx) + # We want subscripting to be read-only so we do not define a __setitem__ method + + def fit_file(self, motion_correct=False, indices=None, include_eval=False, output_dir=None, return_mc=False): """ This method packages the analysis pipeline (motion correction, memory mapping, patch based CNMF processing and component evaluation) in a @@ -307,8 +332,13 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False): perform analysis only on a part of the FOV include_eval (bool) flag for performing component evaluation + output_dir (str) + directory to save the outputs + return_mc (bool) + flag to return motion correction object Returns: cnmf object with the current estimates + (optional) motion correction object """ logger = logging.getLogger("caiman") if indices is None: @@ -321,6 +351,13 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False): logger.error(f"Error: File not found, with file list:\n{fnames[0]}") raise Exception('File not found!') + caiman_temp = os.environ.get("CAIMAN_TEMP") + if output_dir is not None: + # update CAIMAN_TEMP to point to `output_dir` + output_dir = pathlib.Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + os.environ["CAIMAN_TEMP"] = str(output_dir) + base_name = pathlib.Path(fnames[0]).stem + "_memmap_" if extension == '.mmap': fname_new = fnames[0] @@ -363,17 +400,16 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False): fit_cnm = self.fit(images, indices=indices) Cn = caiman.summary_images.local_correlations(images[::max(T//1000, 1)], swap_dim=False) Cn[np.isnan(Cn)] = 0 - fit_cnm.save(fname_new[:-5] + '_init.hdf5') - #fit_cnm.params.change_params({'p': self.params.get('preprocess', 'p')}) - # RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution + fname_init_hdf5 = fname_new[:-5] + '_init.hdf5' + fit_cnm.save(fname_init_hdf5) + # Rerun seeded CNMF on accepted patches to refine and perform deconvolution cnm2 = fit_cnm.refit(images, dview=self.dview) cnm2.estimates.evaluate_components(images, cnm2.params, dview=self.dview) - # update object with selected components - #cnm2.estimates.select_components(use_object=True) # Extract DF/F values cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250) cnm2.estimates.Cn = Cn - cnm2.save(cnm2.mmap_file[:-4] + 'hdf5') + fname_hdf5 = cnm2.mmap_file[:-4] + 'hdf5' + cnm2.save(fname_hdf5) # XXX Why are we stopping the cluster here? What started it? Why remove log files? caiman.cluster.stop_server(dview=self.dview) @@ -381,6 +417,15 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False): for log_file in log_files: os.remove(log_file) + # revert CAIMAN_TEMP to its original value + if caiman_temp is not None: + os.environ["CAIMAN_TEMP"] = caiman_temp + else: + del os.environ["CAIMAN_TEMP"] + + if return_mc & motion_correct: + return cnm2, mc + return cnm2 @@ -403,20 +448,19 @@ def refit(self, images, dview=None): estimates.coordinates = None cnm.estimates = estimates cnm.mmap_file = self.mmap_file - return cnm.fit(images) + cnm.fit(images) + return cnm - def fit(self, images, indices=(slice(None), slice(None))): + def fit(self, images, indices=(slice(None), slice(None))) -> None: """ This method uses the cnmf algorithm to find sources in data. + After it finishes, the C, A, S, b, and f fields will be populated. Args: images : mapped np.ndarray of shape (t,x,y[,z]) containing the images that vary over time. indices: list of slice objects along dimensions (x,y[,z]) for processing only part of the FOV - Returns: - self: updated using the cnmf algorithm with C,A,S,b,f computed according to the given initial values - http://www.cell.com/neuron/fulltext/S0896-6273(15)01084-3 """ @@ -424,11 +468,14 @@ def fit(self, images, indices=(slice(None), slice(None))): # Todo : to compartment if isinstance(indices, slice): indices = [indices] + if isinstance(indices, tuple): indices = list(indices) + indices = [slice(None)] + indices if len(indices) < len(images.shape): indices = indices + [slice(None)]*(len(images.shape) - len(indices)) + dims_orig = images.shape[1:] dims_sliced = images[tuple(indices)].shape[1:] is_sliced = (dims_orig != dims_sliced) @@ -440,7 +487,6 @@ def fit(self, images, indices=(slice(None), slice(None))): T = images.shape[0] self.params.set('online', {'init_batch': T}) self.dims = images.shape[1:] - #self.params.data['dims'] = images.shape[1:] 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): @@ -448,26 +494,18 @@ def fit(self, images, indices=(slice(None), slice(None))): logger.info((T,) + self.dims) - # Make sure filename is pointed correctly (numpy sets it to None sometimes) + # Make sure filename is correctly set (numpy sets it to None sometimes) try: Y.filename = images.filename Yr.filename = images.filename self.mmap_file = images.filename - except AttributeError: # if no memmapping cause working with small data + except AttributeError: # if no memmapping because we're working with small data pass - # update/set all options that depend on data dimensions - # number of rows, columns [and depths] - # self.params.set('spatial', {'medw': (3,) * len(self.dims), - # 'se': np.ones((3,) * len(self.dims), dtype=np.uint8), - # 'ss': np.ones((3,) * len(self.dims), dtype=np.uint8) - # }) - logger.info(f"Using {self.params.get('patch', 'n_processes')} processes") # FIXME The code below is really ugly and it's hard to tell if it's doing the right thing if self.params.get('preprocess', 'n_pixels_per_process') is None: - avail_memory_per_process = psutil.virtual_memory()[ - 1] / 2.**30 / self.params.get('patch', 'n_processes') + avail_memory_per_process = psutil.virtual_memory()[1] / 2.**30 / self.params.get('patch', 'n_processes') mem_per_pix = 3.6977678498329843e-09 npx_per_proc = int(avail_memory_per_process / 8. / mem_per_pix / T) npx_per_proc = int(np.minimum(npx_per_proc, np.prod(self.dims) // self.params.get('patch', 'n_processes'))) @@ -495,29 +533,27 @@ def fit(self, images, indices=(slice(None), slice(None))): if self.remove_very_bad_comps: - logger.info('removing bad components : ') + logger.info('Removing bad components') final_frate = 10 r_values_min = 0.5 # threshold on space consistency fitness_min = -15 # threshold on time variability fitness_delta_min = -15 Npeaks = 10 traces = np.array(self.C) - logger.info('estimating the quality...') + logger.info('Estimating component qualities...') idx_components, idx_components_bad, fitness_raw,\ fitness_delta, r_values = estimate_components_quality( traces, Y, self.estimates.A, self.estimates.C, self.estimates.b, self.estimates.f, final_frate=final_frate, Npeaks=Npeaks, r_values_min=r_values_min, fitness_min=fitness_min, fitness_delta_min=fitness_delta_min, return_all=True, N=5) - logger.info(('Keeping ' + str(len(idx_components)) + - ' and discarding ' + str(len(idx_components_bad)))) + logger.info(f"Keeping {len(idx_components)} components and discarding {len(idx_components_bad)} components") self.estimates.C = self.estimates.C[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() - - return self + return logger.info('update spatial ...') self.update_spatial(Yr, use_init=True) @@ -545,9 +581,6 @@ def fit(self, images, indices=(slice(None), slice(None))): self.params.set('temporal', {'p': self.params.get('preprocess', 'p')}) logger.info('update temporal ...') self.update_temporal(Yr, use_init=False) - # else: - # todo : ask for those.. - # C, f, S, bl, c1, neurons_sn, g1, YrA, lam = self.estimates.C, self.estimates.f, self.estimates.S, self.estimates.bl, self.estimates.c1, self.estimates.neurons_sn, self.estimates.g, self.estimates.YrA, self.estimates.lam # embed in the whole FOV if is_sliced: @@ -576,6 +609,12 @@ def fit(self, images, indices=(slice(None), slice(None))): raise Exception( 'You need to provide a memory mapped file as input if you use patches!!') + # We sometimes need to investigate what changes between before run_CNMF_patches and after that/before we update + # components. This code block here is ready to uncomment for such debugging. + #print("D: About to run run_CNMF_patches(), entering a shell. Will open a shell again afterwards") + #import code + #code.interact(local=dict(globals(), **locals()) ) + self.estimates.A, self.estimates.C, self.estimates.YrA, self.estimates.b, self.estimates.f, \ self.estimates.sn, self.estimates.optional_outputs = run_CNMF_patches( images.filename, self.dims + (T,), self.params, @@ -585,8 +624,14 @@ def fit(self, images, indices=(slice(None), slice(None))): del_duplicates=self.params.get('patch', 'del_duplicates'), indices=indices) + #print("D: Finished with run_CNMF_patches(), self.estimates.* are populated. Next step would be update_temporal() but first: Entering a shell.") + #code.interact(local=dict(globals(), **locals()) ) + + if (self.estimates.b is not None) and (len(list(np.where(~self.estimates.b.any(axis=0))[0])) > 0): # If any of the background ended up completely empty + raise Exception("After run_CNMF_patches(), one or more of the background components is empty. Please restart analysis with init/nb set to a lower value") + self.estimates.bl, self.estimates.c1, self.estimates.g, self.estimates.neurons_sn = None, None, None, None - logger.info("merging") + logger.info("Merging") self.estimates.merged_ROIs = [0] @@ -594,7 +639,7 @@ def fit(self, images, indices=(slice(None), slice(None))): if self.params.get('patch', 'nb_patch') > 0: while len(self.estimates.merged_ROIs) > 0: - self.merge_comps(Yr, mx=np.Inf, fast_merge=True) + self.merge_comps(Yr, mx=np.inf, fast_merge=True) logger.info("update temporal") self.update_temporal(Yr, use_init=False) @@ -607,18 +652,15 @@ def fit(self, images, indices=(slice(None), slice(None))): self.update_temporal(Yr, use_init=False) else: while len(self.estimates.merged_ROIs) > 0: - self.merge_comps(Yr, mx=np.Inf, fast_merge=True) - #if len(self.estimates.merged_ROIs) > 0: - #not_merged = np.setdiff1d(list(range(len(self.estimates.YrA))), - # np.unique(np.concatenate(self.estimates.merged_ROIs))) - #self.estimates.YrA = np.concatenate([self.estimates.YrA[not_merged], - # np.array([self.estimates.YrA[m].mean(0) for ind, m in enumerate(self.estimates.merged_ROIs) if not self.empty_merged[ind]])]) + self.merge_comps(Yr, mx=np.inf, fast_merge=True) + if self.params.get('init', 'nb') == 0: self.estimates.W, self.estimates.b0 = compute_W( Yr, self.estimates.A.toarray(), self.estimates.C, self.dims, self.params.get('init', 'ring_size_factor') * self.params.get('init', 'gSiz')[0], ssub=self.params.get('init', 'ssub_B')) + if len(self.estimates.C): self.deconvolve() self.estimates.C = self.estimates.C.astype(np.float32) @@ -626,14 +668,12 @@ def fit(self, images, indices=(slice(None), slice(None))): self.estimates.S = self.estimates.C else: while len(self.estimates.merged_ROIs) > 0: - self.merge_comps(Yr, mx=np.Inf) + self.merge_comps(Yr, mx=np.inf) - logger.info("update temporal") + logger.info("Updating temporal components") self.update_temporal(Yr, use_init=False) self.estimates.normalize_components() - return self - def save(self, filename): '''save object in hdf5 file format @@ -669,7 +709,7 @@ def remove_components(self, ind_rm): self.params.get('online', 'expected_comps')) self.params.set('online', {'expected_comps': expected_comps}) - def compute_residuals(self, Yr): + def compute_residuals(self, Yr) -> None: """ Compute residual trace for each component (variable YrA). WARNING: At the moment this method is valid only for the 2p processing @@ -704,11 +744,8 @@ def compute_residuals(self, Yr): self.estimates.YrA = (YA - (AA.T.dot(Cf)).T)[:, :self.estimates.A.shape[-1]].T self.estimates.R = self.estimates.YrA - return self - - def deconvolve(self, p=None, method_deconvolution=None, bas_nonneg=None, - noise_method=None, optimize_g=0, s_min=None, **kwargs): + noise_method=None, optimize_g=0, s_min=None, **kwargs) -> None: """Performs deconvolution on already extracted traces using constrained foopsi. """ @@ -744,10 +781,7 @@ def deconvolve(self, p=None, method_deconvolution=None, bas_nonneg=None, else: results = list(map(constrained_foopsi_parallel, args_in)) - if sys.version_info >= (3, 0): - results = list(zip(*results)) - else: # python 2 - results = zip(*results) + results = list(zip(*results)) order = list(results[7]) self.estimates.C = np.stack([results[0][i] for i in order]) @@ -758,106 +792,8 @@ def deconvolve(self, p=None, method_deconvolution=None, bas_nonneg=None, 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 - - def HALS4traces(self, Yr, groups=None, use_groups=False, order=None, - update_bck=True, bck_non_neg=True, **kwargs): - """Solves C, f = argmin_C ||Yr-AC-bf|| using block-coordinate decent. - Can use groups to update non-overlapping components in parallel or a - specified order. - - Args: - Yr : np.array (possibly memory mapped, (x,y,[,z]) x t) - Imaging data reshaped in matrix format - - groups : list of sets - grouped components to be updated simultaneously - - use_groups : bool - flag for using groups - - order : list - Update components in that order (used if nonempty and groups=None) - update_bck : bool - Flag for updating temporal background components - - bck_non_neg : bool - Require temporal background to be non-negative - - Returns: - self (updated values for self.estimates.C, self.estimates.f, self.estimates.YrA) - """ - if update_bck: - Ab = scipy.sparse.hstack([self.estimates.b, self.estimates.A]).tocsc() - try: - Cf = np.vstack([self.estimates.f, self.estimates.C + self.estimates.YrA]) - except(): - Cf = np.vstack([self.estimates.f, self.estimates.C]) - else: - Ab = self.estimates.A - try: - Cf = self.estimates.C + self.estimates.YrA - except(): - Cf = self.estimates.C - Yr = Yr - self.estimates.b.dot(self.estimates.f) - if (groups is None) and use_groups: - 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) # 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')] - else: - self.estimates.f = noisyC[:self.params.get('init', 'nb')] - self.estimates.C = C[self.params.get('init', 'nb'):] - self.estimates.YrA = noisyC[self.params.get('init', 'nb'):] - self.estimates.C - else: - self.estimates.C = C - self.estimates.YrA = noisyC - self.estimates.C - return self - - def HALS4footprints(self, Yr, update_bck=True, num_iter=2): - """Uses hierarchical alternating least squares to update shapes and - background - - Args: - Yr: np.array (possibly memory mapped, (x,y,[,z]) x t) - Imaging data reshaped in matrix format - - update_bck: bool - flag for updating spatial background components - - num_iter: int - number of iterations - - Returns: - self (updated values for self.estimates.A and self.estimates.b) - """ - if update_bck: - Ab = np.hstack([self.estimates.b, self.estimates.A.toarray()]) - try: - Cf = np.vstack([self.estimates.f, self.estimates.C + self.estimates.YrA]) - except(): - Cf = np.vstack([self.estimates.f, self.estimates.C]) - else: - Ab = self.estimates.A.toarray() - try: - Cf = self.estimates.C + self.estimates.YrA - except(): - Cf = self.estimates.C - Yr = Yr - self.estimates.b.dot(self.estimates.f) - 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')] - else: - self.estimates.A = scipy.sparse.csc_matrix(Ab) - - return self - - def update_temporal(self, Y, use_init=True, **kwargs): + def update_temporal(self, Y, use_init=True, **kwargs) -> None: """Updates temporal components Args: @@ -869,24 +805,19 @@ def update_temporal(self, Y, use_init=True, **kwargs): pr = inspect.signature(self.update_temporal) params = [k for k, v in pr.parameters.items() if '=' in str(v)] kw2 = {k: lc[k] for k in params} - try: - kwargs_new = {**kw2, **kwargs} - except(): # python 2.7 - kwargs_new = kw2.copy() - kwargs_new.update(kwargs) + kwargs_new = {**kw2, **kwargs} self.params.set('temporal', kwargs_new) - self.estimates.C, self.estimates.A, self.estimates.b, self.estimates.f, self.estimates.S, \ self.estimates.bl, self.estimates.c1, self.estimates.neurons_sn, \ self.estimates.g, self.estimates.YrA, self.estimates.lam = update_temporal_components( Y, self.estimates.A, self.estimates.b, self.estimates.C, self.estimates.f, dview=self.dview, **self.params.get_group('temporal')) self.estimates.R = self.estimates.YrA - return self - def update_spatial(self, Y, use_init=True, **kwargs): + def update_spatial(self, Y, use_init=True, **kwargs) -> None: """Updates spatial components + modifies values self.estimates.A, self.estimates.b possibly self.estimates.C, self.estimates.f Args: Y: np.array (d1*d2) x T @@ -894,19 +825,12 @@ def update_spatial(self, Y, use_init=True, **kwargs): use_init: bool use Cin, f_in for computing A, b otherwise use C, f - Returns: - self - modified values self.estimates.A, self.estimates.b possibly self.estimates.C, self.estimates.f """ lc = locals() pr = inspect.signature(self.update_spatial) params = [k for k, v in pr.parameters.items() if '=' in str(v)] kw2 = {k: lc[k] for k in params} - try: - kwargs_new = {**kw2, **kwargs} - except(): # python 2.7 - kwargs_new = kw2.copy() - kwargs_new.update(kwargs) + kwargs_new = {**kw2, **kwargs} self.params.set('spatial', kwargs_new) for key in kwargs_new: if hasattr(self, key): @@ -916,9 +840,7 @@ def update_spatial(self, Y, use_init=True, **kwargs): b_in=self.estimates.b, dview=self.dview, sn=self.estimates.sn, dims=self.dims, **self.params.get_group('spatial')) - return self - - def merge_comps(self, Y, mx=50, fast_merge=True): + def merge_comps(self, Y, mx=50, fast_merge=True) -> None: """merges components """ self.estimates.A, self.estimates.C, self.estimates.nr, self.estimates.merged_ROIs, self.estimates.S, \ @@ -931,9 +853,7 @@ def merge_comps(self, Y, mx=50, fast_merge=True): g=self.estimates.g, thr=self.params.get('merging', 'merge_thr'), mx=mx, fast_merge=fast_merge, merge_parallel=self.params.get('merging', 'merge_parallel')) - return self - - def initialize(self, Y, **kwargs): + def initialize(self, Y, **kwargs) -> None: """Component initialization """ self.params.set('init', kwargs) @@ -957,8 +877,6 @@ def initialize(self, Y, **kwargs): self.estimates = estim - return self - def preprocess(self, Yr): """ Examines data to remove corrupted pixels and computes the noise level @@ -969,6 +887,7 @@ def preprocess(self, Yr): 2d array of data (pixels x timesteps) typically in memory mapped form """ + # TODO Weird that this returns Yr Yr, self.estimates.sn, self.estimates.g, self.estimates.psx = preprocess_data( Yr, dview=self.dview, **self.params.get_group('preprocess')) return Yr diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py index d92945a21..cf296071c 100644 --- a/caiman/source_extraction/cnmf/estimates.py +++ b/caiman/source_extraction/cnmf/estimates.py @@ -167,11 +167,35 @@ def __init__(self, A=None, b=None, C=None, f=None, R=None, dims=None): self.A_thr = None self.discarded_components = None - + def __str__(self): + ret = f"Caiman CNMF Estimates Object. subfields:{list(self.__dict__.keys())}" + if hasattr(self, 'A') and self.A is not None: + ret += f" A.shape={self.A.shape}" + if hasattr(self, 'b') and self.b is not None: + ret += f" b.shape={self.b.shape}" + if hasattr(self, 'C') and self.C is not None: + ret += f" C.shape={self.C.shape}" + return ret + + def __repr__(self): + ret = f"Caiman CNMF Estimates Object" + if hasattr(self, 'A') and self.A is not None: + ret += f" A.shape={self.A.shape}" + if hasattr(self, 'b') and self.b is not None and len(self.b.shape) > 1: + ret += f" bg components={self.b.shape[1]}" + if hasattr(self, 'C') and self.C is not None: + ret += f" C.shape={self.C.shape}" + ret += " Use str() for more details" + return ret + + + def __getitem__(self, idx): + return getattr(self, idx) + # We want subscripting to be read-only so we do not define a __setitem__ method def plot_contours(self, img=None, idx=None, thr_method='max', thr=0.2, display_numbers=True, params=None, - cmap='viridis'): + cmap='viridis') -> None: """view contours of all spatial footprints. Args: @@ -226,10 +250,9 @@ def plot_contours(self, img=None, idx=None, thr_method='max', display_numbers=display_numbers, cmap=cmap) plt.title('Rejected Components') - return self def plot_contours_nb(self, img=None, idx=None, thr_method='max', - thr=0.2, params=None, line_color='white', cmap='viridis'): + thr=0.2, params=None, line_color='white', cmap='viridis') -> None: """view contours of all spatial footprints (notebook environment). Args: @@ -309,9 +332,8 @@ def plot_contours_nb(self, img=None, idx=None, thr_method='max', print("Using non-interactive plot as fallback") self.plot_contours(img=img, idx=idx, thr_method=thr_method, thr=thr, params=params, cmap=cmap) - return self - def view_components(self, Yr=None, img=None, idx=None): + def view_components(self, Yr=None, img=None, idx=None) -> None: """view spatial and temporal components interactively Args: @@ -352,10 +374,9 @@ def view_components(self, Yr=None, img=None, idx=None): r_values=None if self.r_values is None else self.r_values[idx], SNR=None if self.SNR_comp is None else self.SNR_comp[idx], cnn_preds=None if np.sum(self.cnn_preds) in (0, None) else self.cnn_preds[idx]) - return self def nb_view_components(self, Yr=None, img=None, idx=None, - denoised_color=None, cmap='jet', thr=0.99): + denoised_color=None, cmap='jet', thr=0.99) -> None: """view spatial and temporal components interactively in a notebook Args: @@ -407,7 +428,6 @@ def nb_view_components(self, Yr=None, img=None, idx=None, r_values=None if self.r_values is None else self.r_values[idx], SNR=None if self.SNR_comp is None else self.SNR_comp[idx], cnn_preds=None if np.sum(self.cnn_preds) in (0, None) else self.cnn_preds[idx]) - return self def hv_view_components(self, Yr=None, img=None, idx=None, denoised_color=None, cmap='viridis'): @@ -463,7 +483,7 @@ def hv_view_components(self, Yr=None, img=None, idx=None, def nb_view_components_3d(self, Yr=None, image_type='mean', dims=None, max_projection=False, axis=0, - denoised_color=None, cmap='jet', thr=0.9): + denoised_color=None, cmap='jet', thr=0.9) -> None: """view spatial and temporal components interactively in a notebook (version for 3d data) @@ -515,8 +535,6 @@ def nb_view_components_3d(self, Yr=None, image_type='mean', dims=None, max_projection=max_projection, axis=axis, thr=thr, denoised_color=denoised_color, cmap=cmap) - return self - def make_color_movie(self, imgs, q_max=99.75, q_min=2, gain_res=1, magnification=1, include_bck=True, frame_range=slice(None, None, None), @@ -739,7 +757,7 @@ def compute_background(self, Yr): (-1, B.shape[-1]), order='F') return B - def compute_residuals(self, Yr): + def compute_residuals(self, Yr) -> None: """compute residual for each component (variable R) Args: @@ -771,11 +789,9 @@ def compute_residuals(self, Yr): AA = Ab.T.dot(Ab) * nA2_inv_mat self.R = (YA - (AA.T.dot(Cf)).T)[:, :self.A.shape[-1]].T - return self - def detrend_df_f(self, quantileMin=8, frames_window=500, flag_auto=True, use_fast=False, use_residuals=True, - detrend_only=False): + detrend_only=False) -> None: """Computes DF/F normalized fluorescence for the extracted traces. See caiman.source.extraction.utilities.detrend_df_f for details @@ -828,9 +844,8 @@ def detrend_df_f(self, quantileMin=8, frames_window=500, frames_window=frames_window, flag_auto=flag_auto, use_fast=use_fast, detrend_only=detrend_only) - return self - def normalize_components(self): + def normalize_components(self) -> None: """ Normalizes components such that spatial components have l_2 norm 1 """ if 'csc_matrix' not in str(type(self.A)): @@ -865,9 +880,8 @@ def normalize_components(self): nB_inv_mat = scipy.sparse.spdiags(1. / (nB + np.finfo(np.float32).eps), 0, nB.shape[0], nB.shape[0]) self.b = self.b * nB_inv_mat self.f = nB_mat * self.f - return self - def select_components(self, idx_components=None, use_object=False, save_discarded_components=True): + def select_components(self, idx_components=None, use_object=False, save_discarded_components=True) -> None: """Keeps only a selected subset of components and removes the rest. The subset can be either user defined with the variable idx_components or read from the estimates object. The flag use_object determines this @@ -946,9 +960,7 @@ def select_components(self, idx_components=None, use_object=False, save_discarde self.idx_components = None self.idx_components_bad = None - return self - - def restore_discarded_components(self): + def restore_discarded_components(self) -> None: ''' Recover components that are filtered out with the select_components method ''' logger = logging.getLogger("caiman") @@ -977,7 +989,7 @@ def restore_discarded_components(self): self.nr = self.A.shape[-1] - def evaluate_components_CNN(self, params, neuron_class=1): + def evaluate_components_CNN(self, params, neuron_class=1) -> None: """Estimates the quality of inferred spatial components using a pretrained CNN classifier. @@ -998,9 +1010,8 @@ class label for neuron shapes predictions = evaluate_components_CNN(self.A, dims, gSig)[0] self.cnn_preds = predictions[:, neuron_class] self.idx_components = np.where(self.cnn_preds >= min_cnn_thr)[0] - return self - def evaluate_components(self, imgs, params, dview=None): + def evaluate_components(self, imgs, params, dview=None) -> None: """Computes the quality metrics for each component and stores the indices of the components that pass user specified thresholds. The various thresholds and parameters can be passed as inputs. If left @@ -1075,9 +1086,8 @@ def evaluate_components(self, imgs, params, dview=None): np.setdiff1d(self.idx_components, idx_ecc)) self.idx_components = np.intersect1d(self.idx_components, idx_ecc) - return self - def filter_components(self, imgs, params, new_dict={}, dview=None, select_mode:str='All'): + def filter_components(self, imgs, params, new_dict={}, dview=None, select_mode:str='All') -> None: """ Filters components based on given thresholds without re-computing the quality metrics. If the quality metrics are not present then it @@ -1165,8 +1175,6 @@ def filter_components(self, imgs, params, new_dict={}, dview=None, select_mode:s self.idx_components_bad = np.array(np.setdiff1d(range(len(self.SNR_comp)),self.idx_components)) - return self - def deconvolve(self, params, dview=None, dff_flag=False): ''' performs deconvolution on the estimated traces using the parameters specified in params. Deconvolution on detrended and normalized (DF/F) diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index f1ff6755c..96a1549d0 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -148,7 +148,8 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter min_corr=0.8, min_pnr=10, seed_method='auto', ring_size_factor=1.5, center_psf=False, ssub_B=2, init_iter=2, remove_baseline = True, SC_kernel='heat', SC_sigma=1, SC_thr=0, SC_normalize=True, SC_use_NN=False, - SC_nnn=20, lambda_gnmf=1, snmf_l1_ratio:float=0.0): + SC_nnn=20, lambda_gnmf=1, snmf_l1_ratio:float=0.0, + greedyroi_nmf_init_method:str='nndsvdar', greedyroi_nmf_max_iter:int=200): """ Initialize components. This function initializes the spatial footprints, temporal components, and background which are then further refined by the CNMF iterations. There are four @@ -168,6 +169,12 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter It is also by default followed by hierarchical alternative least squares (HALS) NMF. Optional use of spatio-temporal downsampling to boost speed. + TODO: Parameters are a mess for this and threading new parameters down to the called code + is awful. In the future we should both clean up the parameters object and just pass it in, + either whole or in some limited (or nested-dict) form, rather than have this many + arguments. This will need to be a complete overhaul because the way params are passed in now + make the structure of CNMFParams closely tied to the function signature. + Args: Y: np.ndarray d1 x d2 [x d3] x T movie, raw data. @@ -259,6 +266,12 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter snmf_l1_ratio: float Used only by sparse NMF, passed to NMF call + greedyroi_nmf_init_method + If greedy_roi is used, this specifies the NMF init method + + greedyroi_nmf_max_iter + If greedy_roi is used, this specifies the max_iter to NMF + Returns: Ain: np.ndarray (d1 * d2 [ * d3]) x K , spatial filter of each neuron. @@ -318,11 +331,11 @@ 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 - logger.info('Roi Initialization...') + logger.info(f'Roi Initialization with method {method}...') 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, seed_method=seed_method) + rolling_sum=rolling_sum, rolling_length=rolling_length, seed_method=seed_method, nmf_overrides={'init_method': greedyroi_nmf_init_method, 'max_iter': greedyroi_nmf_max_iter}) if use_hals: logger.info('Refining Components using HALS NMF iterations') @@ -352,14 +365,13 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter SC_sigma=SC_sigma, SC_use_NN=SC_use_NN, SC_nnn=SC_nnn, SC_normalize=SC_normalize, SC_thr=SC_thr) - elif method == 'pca_ica': + elif method == 'pca_ica': # This code may be untested Ain, Cin, _, b_in, f_in = ICA_PCA( Y_ds, nr=K, sigma_smooth=sigma_smooth_snmf, truncate=2, fun='logcosh', tol=1e-10, max_iter=max_iter_snmf, remove_baseline=True, perc_baseline=perc_baseline_snmf, nb=nb) else: - print(method) - raise Exception("Unsupported initialization method") + raise Exception(f"Unsupported initialization method {method}") K = np.shape(Ain)[-1] @@ -431,6 +443,8 @@ def ICA_PCA(Y_ds, nr, sigma_smooth=(.5, .5, .5), truncate=2, fun='logcosh', nb """ + # Does this work well enough to leave it in the codebase? Unclear if this + # works print("not a function to use in the moment ICA PCA \n") m = scipy.ndimage.gaussian_filter(np.transpose( Y_ds, [2, 0, 1]), sigma=sigma_smooth, mode='nearest', truncate=truncate) @@ -583,30 +597,17 @@ def compressedNMF(Y_ds, nr, r_ov=10, max_iter_snmf=500, T, dims = m.shape[0], m.shape[1:] d = np.prod(dims) yr = np.reshape(m, [T, d], order='F') -# L = randomized_range_finder(yr, nr + r_ov, 3) -# R = randomized_range_finder(yr.T, nr + r_ov, 3) -# Yt = L.T.dot(yr).dot(R) -# c_in, a_in = compressive_nmf(Yt, L, R.T, nr) -# C_in = L.dot(c_in) -# A_in = a_in.dot(R.T) -# A_in = A_in.T -# C_in = C_in.T A, C, USV = nnsvd_init(yr, nr, r_ov=r_ov) W_r = np.random.randn(d, nr + r_ov) W_l = np.random.randn(T, nr + r_ov) US = USV[0]*USV[1] YYt = US.dot(USV[2].dot(USV[2].T)).dot(US.T) -# YYt = yr.dot(yr.T) B = YYt.dot(YYt.dot(US.dot(USV[2].dot(W_r)))) PC, _ = np.linalg.qr(B) B = USV[2].T.dot(US.T.dot(YYt.dot(YYt.dot(W_l)))) PA, _ = np.linalg.qr(B) -# mdl = NMF(n_components=nr, verbose=False, init='nndsvd', tol=1e-10, -# max_iter=1) -# C = mdl.fit_transform(yr).T -# A = mdl.components_.T yrPA = yr.dot(PA) yrPC = PC.T.dot(yr) @@ -690,7 +691,7 @@ def graphNMF(Y_ds, nr, max_iter_snmf=500, lambda_gnmf=1, return A_in, C_in, center, b_in, f_in def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, - rolling_sum=False, rolling_length=100, seed_method='auto'): + rolling_sum=False, rolling_length:int=100, seed_method:str='auto', nmf_overrides:dict = None): """ Greedy initialization of spatial and temporal components using spatial Gaussian filtering @@ -716,7 +717,7 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, nb: int Number of background components - rolling_max: boolean + rolling_sum: boolean Detect new components based on a rolling sum of pixel activity (default: True) rolling_length: int @@ -728,6 +729,17 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, if running as notebook 'semi' and 'manual' require a backend that does not inline figures, e.g. %matplotlib tk + nmf_overrides: dict + This provides a mechanism to override arguments to sklearn.decomposition.NMF() + Right now, the accepted keys are: + 'max_iter' - default number of iterations requested. Defaults to 200 to match the default value in the signature on the sklearn side. + 'init_method' - Default NMF init method. Defaults to 'nndsvdar', the nonnegative double singular value decomposition with + zeroes replaced with small random values. + Some users report better results with 'random', which uses nonnegative random matrices + scaled with sqrt(X.mean() / n_components) + See the docs for sklearn.decomposition.NMF() for other init methods. + Please reach out if you want more options supported. + Returns: A: np.array 2d array of size (# of pixels) x nr with the spatial components. Each column is @@ -739,6 +751,10 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, center: np.array 2d array of size nr x 2 [ or 3] with the components centroids + b_in: (UNDOCUMENTED) + + f_in: (UNDOCUMENTED) + Author: Eftychios A. Pnevmatikakis and Andrea Giovannucci based on a matlab implementation by Yuanjun Gao Simons Foundation, 2015 @@ -749,7 +765,7 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, """ logger = logging.getLogger("caiman") - logger.info("Greedy initialization of spatial and temporal components using spatial Gaussian filtering") + logger.info(f"Greedy initialization of spatial and temporal components using spatial Gaussian filtering. Params={nmf_overrides}") d = np.shape(Y) Y[np.isnan(Y)] = 0 med = np.median(Y, axis=-1) @@ -841,6 +857,7 @@ def greedyROI(Y, nr=30, gSig=[5, 5], gSiz=[11, 11], nIter=5, kernel=None, nb=1, if len(center): plt.scatter(*np.transpose(center)[::-1], c='b') plt.axis('off') + # TODO: Move this graphics code out of the core algorithms; figure out right way to do that plt.suptitle( 'Click to add component. Click again on it to remove it. Press any key to update figure. Add more components, or press any key again when done.') centers = [] @@ -926,8 +943,19 @@ def onclick(event): res = np.reshape(Y, (np.prod(d[0:-1]), d[-1]), order='F') + med.flatten(order='F')[:, None] -# model = NMF(n_components=nb, init='random', random_state=0) - model = NMF(n_components=nb, init='nndsvdar') + + # Defaults + nmf_init_method = 'nndsvdar' + nmf_max_iter = 200 + # Apply overrides + if nmf_overrides is not None and len(nmf_overrides) > 1: + if 'max_iter' in nmf_overrides: + nmf_max_iter = int(nmf_overrides['max_iter']) + if 'init_method' in nmf_overrides: + nmf_init_method = nmf_overrides['init_method'] + + model = NMF(n_components=nb, max_iter=nmf_max_iter, init=nmf_init_method) + b_in = model.fit_transform(np.maximum(res, 0)).astype(np.float32) f_in = model.components_.astype(np.float32) @@ -1253,7 +1281,7 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p logger.info('Merging components') A, C = caiman.source_extraction.cnmf.merging.merge_components( B, A, [], C, None, [], C, [], o, options['spatial'], - dview=None, thr=options['merging']['merge_thr'], mx=np.Inf, fast_merge=True)[:2] + dview=None, thr=options['merging']['merge_thr'], mx=np.inf, fast_merge=True)[:2] A = A.astype(np.float32) C = C.astype(np.float32) logger.info('Updating spatial components') @@ -1302,7 +1330,7 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p logger.info('Merging components') A, C = caiman.source_extraction.cnmf.merging.merge_components( B, A, [], C, None, [], C, [], o, options['spatial'], - dview=None, thr=options['merging']['merge_thr'], mx=np.Inf, fast_merge=True)[:2] + dview=None, thr=options['merging']['merge_thr'], mx=np.inf, fast_merge=True)[:2] A = A.astype(np.float32) C = C.astype(np.float32) logger.info('Updating spatial components') @@ -1333,7 +1361,6 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p 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 = np.linalg.lstsq(b_in, B)[0] else: b_in, s_in, f_in = spr.linalg.svds(B, k=nb) diff --git a/caiman/source_extraction/cnmf/map_reduce.py b/caiman/source_extraction/cnmf/map_reduce.py index 178f3e74e..14153f9d0 100644 --- a/caiman/source_extraction/cnmf/map_reduce.py +++ b/caiman/source_extraction/cnmf/map_reduce.py @@ -109,7 +109,7 @@ def cnmf_patches(args_in): cnm = CNMF(n_processes=1, params=opts) - cnm = cnm.fit(images) + 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, diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index 55b2828ba..6923bf4d2 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/caiman/source_extraction/cnmf/online_cnmf.py @@ -115,6 +115,12 @@ def __init__(self, params=None, estimates=None, path=None, dview=None, Ain=None) self.dview = dview if Ain is not None: self.estimates.A = Ain + if self.params.motion['splits_rig'] > self.params.online['init_batch']/2: + raise Exception("In params, online.init_batch and motion.num_frames_split have incompatible values; consider increasing online.init_batch to be a small multiple of the other") + # See issue #1483; it would actually be better to change initialisation so it ignores splits_rig (either using a default + # value or providing an alternate parameter for that), but that's a much more intrusive change and would potentially + # change things for code/notebooks that've worked for a long time; we should save such changes for a major rewrite + # (if someone takes a particular interest in that). @profile def _prepare_object(self, Yr, T, new_dims=None, idx_components=None): @@ -421,6 +427,7 @@ def fit_next(self, t, frame_in, num_iters_hals=3): num_iters_hals: int, optional maximal number of iterations for HALS (NNLS via blockCD) """ + # FIXME This whole function is overly complex; should rewrite it for legibility logger = logging.getLogger("caiman") t_start = time() @@ -490,7 +497,6 @@ def fit_next(self, t, frame_in, num_iters_hals=3): t_new = time() num_added = 0 if self.params.get('online', 'update_num_comps'): - if self.params.get('online', 'use_corr_img'): corr_img_mode = 'simple' #'exponential' # 'cumulative' self.estimates.corr_img = caiman.summary_images.update_local_correlations( @@ -523,6 +529,7 @@ def fit_next(self, t, frame_in, num_iters_hals=3): else: g_est = 0 use_corr = self.params.get('online', 'use_corr_img') + # FIXME The next statement is really hard to read (self.estimates.Ab, Cf_temp, self.estimates.Yres_buf, self.estimates.rho_buf, self.estimates.CC, self.estimates.CY, self.ind_A, self.estimates.sv, self.estimates.groups, self.estimates.ind_new, self.ind_new_all, @@ -2440,7 +2447,7 @@ def initialize_movie_online(Y, K, gSig, rf, stride, base_name, update_num_comps=True, rval_thr=rval_thr_online, thresh_fitness_delta=thresh_fitness_delta_online, thresh_fitness_raw=thresh_fitness_raw_online, batch_update_suff_stat=True, max_comp_update_shape=5) - cnm_init = cnm_init.fit(images) + cnm_init.fit(images) A_tot = cnm_init.A C_tot = cnm_init.C YrA_tot = cnm_init.YrA @@ -2475,7 +2482,7 @@ def initialize_movie_online(Y, K, gSig, rf, stride, base_name, update_num_comps=True, rval_thr=rval_thr_refine, thresh_fitness_delta=thresh_fitness_delta_refine, thresh_fitness_raw=thresh_fitness_raw_refine, batch_update_suff_stat=True, max_comp_update_shape=5) - cnm_refine = cnm_refine.fit(images) + cnm_refine.fit(images) A, C, b, f, YrA = cnm_refine.A, cnm_refine.C, cnm_refine.b, cnm_refine.f, cnm_refine.YrA diff --git a/caiman/source_extraction/cnmf/params.py b/caiman/source_extraction/cnmf/params.py index 36c9996bb..ff1ae3de0 100644 --- a/caiman/source_extraction/cnmf/params.py +++ b/caiman/source_extraction/cnmf/params.py @@ -215,6 +215,14 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), gSiz: [int, int], default: [int(round((x * 2) + 1)) for x in gSig], half-size of bounding box for each neuron + greedyroi_nmf_init_method: str + When greedyROI is used, this is provided to sklearn's NMF() as the init method. Usually nndsvdar (the default) + is fine, but in some cases random or some other init is preferable; see the sklearn docs for your choices + + greedyroi_nmf_max_iter: int + When greedyROI is used, this is provided to sklearn's NMF() as the max number of iterations; the ideal value + of this partly depends on the init method. See the sklearn docs for guidance. + init_iter: int, default: 2 number of iterations during corr_pnr (1p) initialization @@ -716,6 +724,8 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'gSig': gSig, # size of bounding box 'gSiz': gSiz, + 'greedyroi_nmf_init_method': 'nndsvdar', # init method used in calls to NMF if geedy_roi method for component initialisation is used (offline or online) + 'greedyroi_nmf_max_iter': 200, # max_iter used in calls to NMF if greedy_roi method for component initialisation is used (online or offline) 'init_iter': init_iter, 'kernel': None, # user specified template for greedyROI 'lambda_gnmf': 1, # regularization weight for graph NMF @@ -875,8 +885,8 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'pw_rigid': False, # flag for performing pw-rigid motion correction 'shifts_interpolate': False, # interpolate shifts based on patch locations instead of resizing 'shifts_opencv': True, # flag for applying shifts using cubic interpolation (otherwise FFT) - 'splits_els': 14, # number of splits across time for pw-rigid registration - 'splits_rig': 14, # number of splits across time for rigid registration + 'splits_els': 14, # number of splits across time for pw-rigid registration (usually overridden by code) + 'splits_rig': 14, # number of splits across time for rigid registration (usually overridden by code) 'strides': (96, 96), # how often to start a new patch in pw-rigid registration 'upsample_factor_grid': 4, # motion field upsampling factor during FFT shifts 'use_cuda': False, # flag for using a GPU diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 1ee1ee65a..b0534d75f 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -212,13 +212,14 @@ 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 parallel parrllcomp, len_parrllcomp = caiman.source_extraction.cnmf.utilities.update_order_greedy(AA[:nr, :][:, :nr]) - logger.info("entering the deconvolution ") + logger.info(f"Entering Deconvolution, {C.shape=} {S.shape=}") 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 if np.size(ff) > 0: # Eliminating empty temporal components ff = ff[0] - logger.info(f'removing {len(ff)} empty spatial component(s)') + logger.info(f'removing {len(ff)} empty temporal component(s)') + keep = list(range(A.shape[1])) for i in ff: keep.remove(i) @@ -258,9 +259,6 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, Cin: np.ndarray current estimate of temporal components (K x T) - g: np.ndarray - Global time constant (not used) - bl: np.ndarray baseline for fluorescence trace for each column in A @@ -314,18 +312,21 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, bl: float same as input + YrA: np.ndarray + matrix of spatial component filtered raw data, after all contributions have been removed. + YrA corresponds to the residual trace for each component and is used for faster plotting (K x T) + c1: float same as input - g: float + sn: float same as input - sn: float + g: float same as input - YrA: np.ndarray - matrix of spatial component filtered raw data, after all contributions have been removed. - YrA corresponds to the residual trace for each component and is used for faster plotting (K x T) + lam: (undocumented) + """ logger = logging.getLogger("caiman") @@ -381,7 +382,7 @@ def update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr, f" out of total {nr} temporal components updated") for ii in np.arange(nr, nr + nb): - cc = np.maximum(YrA[:, ii] + Cin[ii], -np.Inf) + cc = np.maximum(YrA[:, ii] + Cin[ii], -np.inf) YrA -= AA[ii, :].T.dot((cc - Cin[ii])[None, :]).T C[ii, :] = cc diff --git a/caiman/source_extraction/volpy/volpy_gui.py b/caiman/source_extraction/volpy/volpy_gui.py old mode 100644 new mode 100755 index ae333dfc1..d8710de38 --- a/caiman/source_extraction/volpy/volpy_gui.py +++ b/caiman/source_extraction/volpy/volpy_gui.py @@ -6,6 +6,13 @@ and spike extraction step of VolPy. @author: @caichangjia """ + +# These imports apparently must come before importing pyqtgraph on some platforms +import PySide6 +from PySide6 import QtWidgets +from PySide6.QtWidgets import QApplication +from PySide6.QtGui import QShortcut + import cv2 import h5py import numpy as np @@ -15,18 +22,15 @@ from pyqtgraph import FileDialog from pyqtgraph.Qt import QtGui from pyqtgraph.parametertree import Parameter, ParameterTree -import PyQt5 -from PyQt5 import QtWidgets -from PyQt5.QtWidgets import QApplication, QShortcut import random from skimage.draw import polygon import sys import caiman as cm from caiman.external.cell_magic_wand import cell_magic_wand_single_point - + os.environ["QT_QPA_PLATFORM_PLUGIN_PATH"] = os.fspath( - Path(PyQt5.__file__).resolve().parent / "Qt5" / "plugins") + Path(PySide6.__file__).resolve().parent / "Qt6" / "plugins") def mouseClickEvent(event): @@ -70,7 +74,7 @@ def mouseClickEvent(event): roi = roi * 1. except: pass - + elif mode == 'CHOOSE NEURONS': pos = img.mapFromScene(event.pos()) p1.clear() @@ -104,7 +108,7 @@ def add(): roi = np.zeros((dims[1], dims[0])) ff = np.array(polygon(np.array(pts)[:,0], np.array(pts)[:,1])) roi[ff[1],[ff[0]]] = 1 - + if len(pts) > 2 : flag = True while flag: @@ -187,7 +191,7 @@ def load_rois(): if (l_ROIs.shape[2], l_ROIs.shape[1]) != dims: print(dims);print(l_ROIs.shape[1:]) raise ValueError('Dimentions of movie and rois do not accord') - + for roi in l_ROIs: flag = True while flag: @@ -213,7 +217,7 @@ def save(): print(ffll[0]) save_ROIs = np.array(list(all_ROIs.values())).copy() save_ROIs = np.flip(save_ROIs, axis=1) - + if os.path.splitext(ffll[0])[1] == '.hdf5': cm.movie(save_ROIs).save(ffll[0]) summary_images = summary_images.transpose([0, 2, 1]) @@ -293,11 +297,13 @@ def overlay(all_ROIs): if __name__ == "__main__": ## Always start by initializing Qt (only once per application) + if sys.platform == 'darwin': + PySide6.QtWidgets.QApplication.setStyle("fusion") app = QApplication(sys.argv) - + ## Define a top-level widget to hold everything w = QtWidgets.QWidget() - + ## Create some widgets to be placed inside hist = pg.HistogramLUTItem() # Contrast/color control win = pg.GraphicsLayoutWidget() @@ -307,11 +313,11 @@ def overlay(all_ROIs): p1 = pg.PlotWidget() neuron_action = ParameterTree() neuron_list = QtWidgets.QListWidget() - + ## Create a grid layout to manage the widgets size and position layout = QtWidgets.QGridLayout() w.setLayout(layout) - + ## Add widgets to the layout in their proper positions layout.addWidget(win, 0, 1) layout.addWidget(p1, 0, 2) @@ -320,27 +326,27 @@ def overlay(all_ROIs): img = pg.ImageItem() p1.addItem(img) hist.setImageItem(img) - + # Add actions - params_action = [{'name': 'LOAD DATA', 'type':'action'}, - {'name': 'LOAD ROIS', 'type':'action'}, - {'name': 'SAVE', 'type':'action'}, - {'name': 'ADD', 'type': 'action'}, - {'name': 'REMOVE', 'type': 'action'}, - {'name': 'SHOW ALL', 'type': 'action'}, - {'name': 'CLEAR', 'type': 'action'}, - {'name': 'IMAGES', 'type': 'list', 'values': ['MEAN','CORR']}, - {'name': 'DISPLAY', 'type': 'list', 'values': ['CONTOUR','SPATIAL FOOTPRINTS']}, - {'name': 'MODE', 'type': 'list', 'values': ['POLYGON','CELL MAGIC WAND', 'CHOOSE NEURONS']}, + params_action = [{'name': 'LOAD DATA', 'type': 'action'}, + {'name': 'LOAD ROIS', 'type': 'action'}, + {'name': 'SAVE', 'type': 'action'}, + {'name': 'ADD', 'type': 'action'}, + {'name': 'REMOVE', 'type': 'action'}, + {'name': 'SHOW ALL', 'type': 'action'}, + {'name': 'CLEAR', 'type': 'action'}, + {'name': 'IMAGES', 'type': 'list', 'limits': ['MEAN','CORR']}, + {'name': 'DISPLAY', 'type': 'list', 'limits': ['CONTOUR', 'SPATIAL FOOTPRINTS']}, + {'name': 'MODE', 'type': 'list', 'limits': ['POLYGON', 'CELL MAGIC WAND', 'CHOOSE NEURONS']}, {'name': 'MAGIC WAND PARAMS', 'type': 'group', 'children': [{'name': 'MIN RADIUS', 'type': 'int', 'value': 4}, {'name': 'MAX RADIUS', 'type': 'int', 'value': 10}, {'name': 'ROUGHNESS', 'type': 'int', 'value': 1}]}] - + pars_action = Parameter.create(name='params_action', type='group', children=params_action) neuron_action.setParameters(pars_action, showTop=False) neuron_action.setWindowTitle('Parameter Action') mode = pars_action.getValues()['MODE'][0] - + # Add event p1.mousePressEvent = mouseClickEvent p1.mouseReleaseEvent = release @@ -358,7 +364,7 @@ def overlay(all_ROIs): shortcut_up = QShortcut(QtGui.QKeySequence("up"), w) shortcut_up.activated.connect(up) neuron_list.itemClicked.connect(show_neuron) - + # Create dictionary for saving F = FileDialog() all_pts = {} @@ -366,12 +372,9 @@ def overlay(all_ROIs): all_ROIs = {} pts = [] pen = pg.mkPen(color=(255, 255, 0), width=4)#, style=Qt.DashDotLine) - + ## Display the widget as a new window w.show() - + ## Start the Qt event loop - app.exec_() - - - + app.exec() diff --git a/caiman/tests/comparison/create_gt.py b/caiman/tests/comparison/create_gt.py index 1954f7c4f..c719cf62e 100644 --- a/caiman/tests/comparison/create_gt.py +++ b/caiman/tests/comparison/create_gt.py @@ -239,7 +239,7 @@ def create(): method_deconvolution='oasis') comp.cnmpatch = copy.copy(cnm) comp.cnmpatch.estimates = None - cnm = cnm.fit(images) + cnm.fit(images) A_tot = cnm.estimates.A C_tot = cnm.estimates.C YrA_tot = cnm.estimates.YrA @@ -286,7 +286,7 @@ def create(): rf=None, stride=None, method_deconvolution='oasis') - cnm = cnm.fit(images) + cnm.fit(images) # DISCARDING A, C, b, f, YrA, sn = cnm.estimates.A, cnm.estimates.C, cnm.estimates.b, cnm.estimates.f, cnm.estimates.YrA, cnm.estimates.sn final_frate = params_movie['final_frate'] diff --git a/caiman/tests/comparison_general.py b/caiman/tests/comparison_general.py index 5361dad28..e7698a2c4 100644 --- a/caiman/tests/comparison_general.py +++ b/caiman/tests/comparison_general.py @@ -242,7 +242,7 @@ def test_general(): method_deconvolution='oasis') comp.cnmpatch = copy.copy(cnm) comp.cnmpatch.estimates = None - cnm = cnm.fit(images) + cnm.fit(images) A_tot = cnm.estimates.A C_tot = cnm.estimates.C YrA_tot = cnm.estimates.YrA @@ -289,7 +289,7 @@ def test_general(): rf=None, stride=None, method_deconvolution='oasis') - cnm = cnm.fit(images) + cnm.fit(images) # DISCARDING A, C, b, f, YrA, sn = cnm.estimates.A, cnm.estimates.C, cnm.estimates.b, cnm.estimates.f, cnm.estimates.YrA, cnm.estimates.sn final_frate = params_movie['final_frate'] diff --git a/caiman/tests/comparison_humans.py b/caiman/tests/comparison_humans.py index 38fcca106..2a3497cd2 100644 --- a/caiman/tests/comparison_humans.py +++ b/caiman/tests/comparison_humans.py @@ -369,7 +369,7 @@ t1 = time.time() print('Starting CNMF') cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) - cnm = cnm.fit(images) + cnm.fit(images) t_patch = time.time() - t1 # %% try: diff --git a/caiman/tests/test_demo.py b/caiman/tests/test_demo.py index 190180012..a1ae9e520 100644 --- a/caiman/tests/test_demo.py +++ b/caiman/tests/test_demo.py @@ -34,7 +34,7 @@ def demo(parallel=False): method_deconvolution='oasis') # FIT images = np.reshape(Yr.T, [T] + list(dims), order='F') - cnm = cnm.fit(images) + cnm.fit(images) if parallel: caiman.cluster.stop_server(dview=dview) diff --git a/caiman/tests/test_sbx.py b/caiman/tests/test_sbx.py index d8609edaa..cba2abf63 100644 --- a/caiman/tests/test_sbx.py +++ b/caiman/tests/test_sbx.py @@ -29,6 +29,7 @@ def test_load_2d(): assert data_2d.ndim == 3, 'Loaded 2D data has wrong dimensionality' assert data_2d.shape == SHAPE_2D, 'Loaded 2D data has wrong shape' assert data_2d.shape == (meta_2d['num_frames'], *meta_2d['frame_size']), 'Shape in metadata does not match loaded data' + assert meta_2d['frame_rate'] == 15.625, 'Frame rate in metadata is incorrect (unidirectional)' npt.assert_array_equal(data_2d[0, 0, :10], [712, 931, 1048, 825, 1383, 882, 601, 798, 1022, 966], 'Loaded 2D data has wrong values') data_2d_movie = cm.load(file_2d) @@ -37,6 +38,23 @@ def test_load_2d(): npt.assert_array_almost_equal(data_2d_movie, data_2d, err_msg='Movie loaded with cm.load has wrong values') +def test_load_2d_bidi(): + file_2d_bidi = os.path.join(TESTDATA_PATH, '2d_sbx_bidi.sbx') + data_2d_bidi = sbx_utils.sbxread(file_2d_bidi) + meta_2d_bidi = sbx_utils.sbx_meta_data(file_2d_bidi) + + assert data_2d_bidi.ndim == 3, 'Loaded 2D bidirectional data has wrong dimensionality' + assert data_2d_bidi.shape == SHAPE_2D, 'Loaded 2D bidirectional data has wrong shape' + assert data_2d_bidi.shape == (meta_2d_bidi['num_frames'], *meta_2d_bidi['frame_size']), 'Shape in metadata does not match loaded data' + assert meta_2d_bidi['frame_rate'] == 31.25, 'Frame rate in metadata is incorrect (bidirectional)' + npt.assert_array_equal(data_2d_bidi[0, 0, :10], [2833, 1538, 1741, 1837, 2079, 2038, 1946, 1631, 2260, 2073], 'Loaded 2D bidirectional data has wrong values') + + data_2d_bidi_movie = cm.load(file_2d_bidi) + assert data_2d_bidi_movie.ndim == data_2d_bidi.ndim, 'Movie loaded with cm.load has wrong dimensionality' + assert data_2d_bidi_movie.shape == data_2d_bidi.shape, 'Movie loaded with cm.load has wrong shape' + npt.assert_array_almost_equal(data_2d_bidi_movie, data_2d_bidi, err_msg='Movie loaded with cm.load has wrong values') + + def test_load_3d(): file_3d = os.path.join(TESTDATA_PATH, '3d_sbx_1.sbx') data_3d = sbx_utils.sbxread(file_3d) @@ -45,6 +63,7 @@ def test_load_3d(): assert data_3d.ndim == 4, 'Loaded 3D data has wrong dimensionality' assert data_3d.shape == SHAPE_3D, 'Loaded 3D data has wrong shape' assert data_3d.shape == (meta_3d['num_frames'], *meta_3d['frame_size'], meta_3d['num_planes']), 'Shape in metadata does not match loaded data' + assert meta_3d['frame_rate'] == 15.625 / meta_3d['num_planes'], 'Frame rate in metadata is incorrect (bidirectional 3D)' npt.assert_array_equal(data_3d[0, 0, :10, 0], [2167, 2525, 1713, 1747, 1887, 1741, 1873, 1244, 1747, 1637], 'Loaded 2D data has wrong values') data_3d_movie = cm.load(file_3d, is3D=True) diff --git a/caiman/tests/test_toydata.py b/caiman/tests/test_toydata.py index c796eaddd..74a6d5adb 100644 --- a/caiman/tests/test_toydata.py +++ b/caiman/tests/test_toydata.py @@ -118,7 +118,7 @@ def pipeline(D, params_dict, name): images = np.reshape(images.T, (T,) + dims, order='F') params.change_params({'data': {'fnames': [mmap_name]}}) - cnm = cnm.fit(images) + cnm.fit(images) # VERIFY HIGH CORRELATION WITH GROUND TRUTH sorting = [np.argmax([np.corrcoef(tc, c)[0, 1] for tc in trueC]) for c in cnm.estimates.C] diff --git a/caiman/utils/sbx_utils.py b/caiman/utils/sbx_utils.py index 48bbba484..610b1be93 100644 --- a/caiman/utils/sbx_utils.py +++ b/caiman/utils/sbx_utils.py @@ -9,13 +9,13 @@ import os import scipy import tifffile -from typing import Iterable, Union, Optional +from typing import Any, Sequence, Union, Optional, cast -DimSubindices = Union[Iterable[int], slice] -FileSubindices = Union[DimSubindices, Iterable[DimSubindices]] # can have inds for just frames or also for y, x, z -ChainSubindices = Union[FileSubindices, Iterable[FileSubindices]] # one to apply to each file, or separate for each file +DimSubindices = Union[Sequence[int], slice] +FileSubindices = Union[DimSubindices, Sequence[DimSubindices]] # can have inds for just frames or also for y, x, z +ChainSubindices = Union[FileSubindices, Sequence[FileSubindices]] # one to apply to each file, or separate for each file -def loadmat_sbx(filename: str) -> dict: +def loadmat_sbx(filename: str) -> dict[str, Any]: """ this wrapper should be called instead of directly calling spio.loadmat @@ -25,7 +25,7 @@ def loadmat_sbx(filename: str) -> dict: """ data_ = scipy.io.loadmat(filename, struct_as_record=False, squeeze_me=True) _check_keys(data_) - return data_ + return data_['info'] def _check_keys(checkdict: dict) -> None: @@ -135,7 +135,7 @@ def sbx_chain_to_tif(filenames: list[str], fileout: str, subindices: Optional[Ch fileout: str filename to save, including the .tif suffix - subindices: Iterable[int] | slice | Iterable[Iterable[int] | slice | tuple[Iterable[int] | slice, ...]] + subindices: Sequence[int] | slice | Sequence[Sequence[int] | slice | tuple[Sequence[int] | slice, ...]] see subindices for sbx_to_tif can specify separate subindices for each file if nested 2 levels deep; X, Y, and Z sizes must match for all files after indexing. @@ -146,21 +146,25 @@ def sbx_chain_to_tif(filenames: list[str], fileout: str, subindices: Optional[Ch subindices = slice(None) # Validate aggressively to avoid failing after waiting to copy a lot of data - if isinstance(subindices, slice) or np.isscalar(subindices[0]): - # One set of subindices to repeat for each file - subindices = [(subindices,) for _ in filenames] + if isinstance(subindices, slice) or isinstance(subindices[0], int) or np.isscalar(subindices[0]): + # One set of subindices over time to repeat for each file + _subindices = [(cast(DimSubindices, subindices),) for _ in filenames] - elif isinstance(subindices[0], slice) or np.isscalar(subindices[0][0]): + elif isinstance(subindices[0], slice) or isinstance(subindices[0][0], int) or np.isscalar(subindices[0][0]): # Interpret this as being an iterable over dimensions to repeat for each file - subindices = [subindices for _ in filenames] + _subindices = [cast(FileSubindices, subindices) for _ in filenames] elif len(subindices) != len(filenames): # Must be a separate subindices for each file; must match number of files - raise Exception('Length of subindices does not match length of file list') + raise Exception('Length of subindices does not match length of file list') + + else: + _subindices = cast(Sequence[FileSubindices], subindices) + del subindices # ensure _subindices replaces subindices from here # Get the total size of the file all_shapes = [sbx_shape(file) for file in filenames] - all_shapes_out = np.stack([_get_output_shape(file, subind)[0] for (file, subind) in zip(filenames, subindices)]) + all_shapes_out = np.stack([_get_output_shape(file, subind)[0] for (file, subind) in zip(filenames, _subindices)]) # Check that X, Y, and Z are consistent for dimname, shapes in zip(('Y', 'X', 'Z'), all_shapes_out.T[1:]): @@ -199,14 +203,15 @@ def sbx_chain_to_tif(filenames: list[str], fileout: str, subindices: Optional[Ch # Now convert each file tif_memmap = tifffile.memmap(fileout, series=0) offset = 0 - for filename, subind, file_N in zip(filenames, subindices, all_n_frames_out): - _sbxread_helper(filename, subindices=subind, channel=channel, out=tif_memmap[offset:offset+file_N], plane=plane, chunk_size=chunk_size) + for filename, subind, file_N in zip(filenames, _subindices, all_n_frames_out): + this_memmap = cast(np.memmap, tif_memmap[offset:offset+file_N]) + _sbxread_helper(filename, subindices=subind, channel=channel, out=this_memmap, plane=plane, chunk_size=chunk_size) offset += file_N del tif_memmap # important to make sure file is closed (on Windows) -def sbx_shape(filename: str, info: Optional[dict] = None) -> tuple[int, int, int, int, int]: +def sbx_shape(filename: str, info: Optional[dict[str, Any]] = None) -> tuple[int, int, int, int, int]: """ Args: filename: str @@ -223,55 +228,41 @@ def sbx_shape(filename: str, info: Optional[dict] = None) -> tuple[int, int, int # Load info if info is None: - info = loadmat_sbx(filename + '.mat')['info'] + info = loadmat_sbx(filename + '.mat') # Image size if 'sz' not in info: info['sz'] = np.array([512, 796]) - - # Scan mode (0 indicates bidirectional) - if 'scanmode' in info and info['scanmode'] == 0: - info['recordsPerBuffer'] *= 2 # Fold lines (multiple subframes per scan) - basically means the frames are smaller and # there are more of them than is reflected in the info file if 'fold_lines' in info and info['fold_lines'] > 0: - if info['recordsPerBuffer'] % info['fold_lines'] != 0: + if info['sz'][0] % info['fold_lines'] != 0: raise Exception('Non-integer folds per frame not supported') - n_folds = round(info['recordsPerBuffer'] / info['fold_lines']) - info['recordsPerBuffer'] = info['fold_lines'] + info['sz'][0] = info['fold_lines'] if 'bytesPerBuffer' in info: + n_folds = round(info['sz'][0] / info['fold_lines']) info['bytesPerBuffer'] /= n_folds - else: - n_folds = 1 + # Defining number of channels/size factor if 'chan' in info: info['nChan'] = info['chan']['nchan'] - factor = 1 # should not be used + elif info['channels'] == 1: + info['nChan'] = 2 else: - if info['channels'] == 1: - info['nChan'] = 2 - factor = 1 - elif info['channels'] == 2: - info['nChan'] = 1 - factor = 2 - elif info['channels'] == 3: - info['nChan'] = 1 - factor = 2 + info['nChan'] = 1 # Determine number of frames in whole file filesize = os.path.getsize(filename + '.sbx') if 'scanbox_version' in info: - if info['scanbox_version'] == 2: - info['max_idx'] = filesize / info['recordsPerBuffer'] / info['sz'][1] * factor / 4 - 1 - elif info['scanbox_version'] == 3: + if info['scanbox_version'] in [2, 3]: info['max_idx'] = filesize / np.prod(info['sz']) / info['nChan'] / 2 - 1 else: raise Exception('Invalid Scanbox version') else: - info['max_idx'] = filesize / info['bytesPerBuffer'] * factor - 1 + info['max_idx'] = filesize / info['bytesPerBuffer'] * (2 // info['nChan']) - 1 n_frames = info['max_idx'] + 1 # Last frame @@ -284,7 +275,7 @@ def sbx_shape(filename: str, info: Optional[dict] = None) -> tuple[int, int, int n_planes = 1 n_frames //= n_planes - x = (int(info['nChan']), int(info['sz'][1]), int(info['recordsPerBuffer']), int(n_planes), int(n_frames)) + x = (int(info['nChan']), int(info['sz'][1]), int(info['sz'][0]), int(n_planes), int(n_frames)) return x @@ -302,7 +293,7 @@ def sbx_meta_data(filename: str): if ext == '.sbx': filename = basename - info = loadmat_sbx(filename + '.mat')['info'] + info = loadmat_sbx(filename + '.mat') meta_data = dict() n_chan, n_x, n_y, n_planes, n_frames = sbx_shape(filename, info) @@ -400,13 +391,14 @@ def _sbxread_helper(filename: str, subindices: FileSubindices = slice(None), cha filename = basename # Normalize so subindices is a list over dimensions - if isinstance(subindices, slice) or np.isscalar(subindices[0]): - subindices = [subindices] + if isinstance(subindices, slice) or isinstance(subindices[0], int) or np.isscalar(subindices[0]): + _subindices = [cast(DimSubindices, subindices)] else: - subindices = list(subindices) + _subindices = list(cast(Sequence[DimSubindices], subindices)) + del subindices # ensure _subindices replaces subindices from here # Load info - info = loadmat_sbx(filename + '.mat')['info'] + info = loadmat_sbx(filename + '.mat') # Get shape (and update info) data_shape = sbx_shape(filename, info) # (chans, X, Y, Z, frames) @@ -414,7 +406,7 @@ def _sbxread_helper(filename: str, subindices: FileSubindices = slice(None), cha is3D = n_planes > 1 # Fill in missing dimensions in subindices - subindices += [slice(None) for _ in range(max(0, 3 + is3D - len(subindices)))] + _subindices += [slice(None) for _ in range(max(0, 3 + is3D - len(_subindices)))] if channel is None: if n_chans > 1: @@ -430,7 +422,7 @@ def _sbxread_helper(filename: str, subindices: FileSubindices = slice(None), cha if frame_size <= 0: raise Exception('Invalid scanbox metadata') - save_shape, subindices = _get_output_shape(data_shape, subindices) + save_shape, _subindices = _get_output_shape(data_shape, _subindices) n_frames_out = save_shape[0] if plane is not None: if len(save_shape) < 4: @@ -447,16 +439,18 @@ def _sbxread_helper(filename: str, subindices: FileSubindices = slice(None), cha if not is3D: # squeeze out singleton plane dim sbx_mmap = sbx_mmap[..., 0] elif plane is not None: # select plane relative to subindices - sbx_mmap = sbx_mmap[..., subindices[-1][plane]] - subindices = subindices[:-1] - inds = np.ix_(*subindices) + sbx_mmap = sbx_mmap[..., _subindices[-1][plane]] + _subindices = _subindices[:-1] + inds = np.ix_(*_subindices) + out_arr: Optional[np.ndarray] = out # widen type + del out # ensure out_arr replaces out from here if chunk_size is None: # load a contiguous block all at once chunk_size = n_frames_out - elif out is None: + elif out_arr is None: # Pre-allocate destination when loading in chunks - out = np.empty(save_shape, dtype=np.uint16) + out_arr = np.empty(save_shape, dtype=np.uint16) n_remaining = n_frames_out offset = 0 @@ -468,23 +462,26 @@ def _sbxread_helper(filename: str, subindices: FileSubindices = slice(None), cha # Note: SBX files store the values strangely, it's necessary to invert each uint16 value to get the correct ones np.invert(chunk, out=chunk) # avoid copying, may be large - if out is None: - out = chunk # avoid copying when loading all data + if out_arr is None: + out_arr = chunk # avoid copying when loading all data else: - out[offset:offset+this_chunk_size] = chunk + out_arr[offset:offset+this_chunk_size] = chunk n_remaining -= this_chunk_size offset += this_chunk_size + if out_arr is None: + raise RuntimeError('Nothing loaded - no frames selected?') + del sbx_mmap # Important to close file (on Windows) - if isinstance(out, np.memmap): - out.flush() - return out + if isinstance(out_arr, np.memmap): + out_arr.flush() + return out_arr -def _interpret_subindices(subindices: DimSubindices, dim_extent: int) -> tuple[Iterable[int], int]: +def _interpret_subindices(subindices: DimSubindices, dim_extent: int) -> tuple[Sequence[int], int]: """ - Given the extent of a dimension in the corresponding recording, obtain an iterable over subindices + Given the extent of a dimension in the corresponding recording, obtain a sequence over subindices and the step size (or 0 if the step size is not uniform). """ logger = logging.getLogger("caiman") @@ -507,15 +504,18 @@ def _interpret_subindices(subindices: DimSubindices, dim_extent: int) -> tuple[I def _get_output_shape(filename_or_shape: Union[str, tuple[int, ...]], subindices: FileSubindices - ) -> tuple[tuple[int, ...], FileSubindices]: + ) -> tuple[tuple[int, ...], tuple[Sequence[int], ...]]: """ Helper to determine what shape will be loaded/saved given subindices Also returns back the subindices with slices transformed to ranges, for convenience """ if isinstance(subindices, slice) or np.isscalar(subindices[0]): - subindices = (subindices,) + _subindices = (cast(DimSubindices, subindices),) + else: + _subindices = cast(Sequence[DimSubindices], subindices) + del subindices # ensure _subindices replaces subindices from here - n_inds = len(subindices) # number of dimensions that are indexed + n_inds = len(_subindices) # number of dimensions that are indexed if isinstance(filename_or_shape, str): data_shape = sbx_shape(filename_or_shape) @@ -529,7 +529,7 @@ def _get_output_shape(filename_or_shape: Union[str, tuple[int, ...]], subindices shape_out = [n_frames, n_y, n_x, n_planes] if is3D else [n_frames, n_y, n_x] subinds_out = [] - for i, (dim, subind) in enumerate(zip(shape_out, subindices)): + for i, (dim, subind) in enumerate(zip(shape_out, _subindices)): iterable_elements = _interpret_subindices(subind, dim)[0] shape_out[i] = len(iterable_elements) subinds_out.append(iterable_elements) diff --git a/caiman/utils/utils.py b/caiman/utils/utils.py index 112f616bf..d94a8eaf3 100644 --- a/caiman/utils/utils.py +++ b/caiman/utils/utils.py @@ -682,7 +682,7 @@ def get_caiman_version() -> tuple[str, str]: class caitimer(contextlib.ContextDecorator): """ This is a simple context manager that you can use like this to get timing information on functions you call: with caiman.utils.utils.caitimer("CNMF fit"): - cnm = cnm.fit(images) + cnm.fit(images) When the context exits it will say how long it was open. Useful for easy function benchmarking """ diff --git a/demos/general/demo_caiman_basic.py b/demos/general/demo_caiman_basic.py index 4731a7a18..be3381ca5 100755 --- a/demos/general/demo_caiman_basic.py +++ b/demos/general/demo_caiman_basic.py @@ -64,7 +64,7 @@ def main(): # Run CaImAn Batch (CNMF) cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) - cnm = cnm.fit_file() + cnm.fit_file() # plot contour plots of components Cns = local_correlations_movie_offline(fnames[0], diff --git a/demos/general/demo_pipeline.py b/demos/general/demo_pipeline.py index 57ef191c4..619703ba2 100755 --- a/demos/general/demo_pipeline.py +++ b/demos/general/demo_pipeline.py @@ -118,7 +118,7 @@ def main(): #opts.change_params({'p': 0}) cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) - cnm = cnm.fit(images) + cnm.fit(images) # plot contours of found components Cns = local_correlations_movie_offline(mc.mmap_file[0], diff --git a/demos/general/demo_pipeline_NWB.py b/demos/general/demo_pipeline_NWB.py index f2983d2fb..79d13b3a5 100755 --- a/demos/general/demo_pipeline_NWB.py +++ b/demos/general/demo_pipeline_NWB.py @@ -165,7 +165,7 @@ def main(): saved_p = opts.preprocess['p'] # Save the deconvolution parameter for later restoration opts.change_params({'preprocess': {'p': 0}, 'temporal': {'p': 0}}) cnm = cnmf.CNMF(n_processes, params=opts, dview=dview) - cnm = cnm.fit(images) + cnm.fit(images) # ALTERNATE WAY TO RUN THE PIPELINE AT ONCE (optional) diff --git a/demos/general/demo_pipeline_voltage_imaging.py b/demos/general/demo_pipeline_voltage_imaging.py index 6ac39faed..4734de079 100755 --- a/demos/general/demo_pipeline_voltage_imaging.py +++ b/demos/general/demo_pipeline_voltage_imaging.py @@ -52,8 +52,8 @@ def main(): # Figure out what data we're working on if cfg.input is None: # If no input is specified, use sample data, downloading if necessary - fnames = [download_demo('demo_voltage_imaging.hdf5', 'volpy')] # XXX do we need to use a separate directory? - path_ROIs = [download_demo('demo_voltage_imaging_ROIs.hdf5', 'volpy')] + fnames = [download_demo('demo_voltage_imaging.hdf5')] # XXX do we need to use a separate directory? + path_ROIs = [download_demo('demo_voltage_imaging_ROIs.hdf5')] else: fnames = cfg.input path_ROIs = cfg.pathinput @@ -74,7 +74,7 @@ def main(): max_deviation_rigid = 3 # maximum deviation allowed for patch with respect to rigid shifts border_nan = 'copy' - params_dict = {'fnames': fnames, + opts_dict = {'fnames': fnames, 'fr': fr, 'pw_rigid': pw_rigid, 'max_shifts': max_shifts, @@ -84,7 +84,7 @@ def main(): 'max_deviation_rigid': max_deviation_rigid, 'border_nan': border_nan} - opts = volparams(params_dict=params_dict) + opts = volparams(params_dict=opts_dict) # play the movie (optional) # playing the movie using opencv. It requires loading the movie in memory. @@ -167,7 +167,7 @@ def main(): elif cfg.method == 'gui_annotation': # run volpy_gui.py file in the caiman/source_extraction/volpy folder - gui_ROIs = caiman_datadir() + '/example_movies/volpy/gui_roi.hdf5' + gui_ROIs = os.path.join(caiman_datadir(), 'example_movies', 'gui_roi.hdf5') with h5py.File(gui_ROIs, 'r') as fl: ROIs = fl['mov'][()] @@ -221,7 +221,7 @@ def main(): 'weight_update': weight_update, 'n_iter': n_iter} - opts.change_params(params_dict=params_dict); + opts.change_params(params_dict=opts_dict); # TRACE DENOISING AND SPIKE DETECTION vpy = VOLPY(n_processes=n_processes, dview=dview, params=opts) diff --git a/demos/notebooks/demo_caiman_cnmf_3D.ipynb b/demos/notebooks/demo_caiman_cnmf_3D.ipynb index 607758214..dd5c64b76 100644 --- a/demos/notebooks/demo_caiman_cnmf_3D.ipynb +++ b/demos/notebooks/demo_caiman_cnmf_3D.ipynb @@ -399,7 +399,7 @@ "source": [ "%%capture\n", "# FIT\n", - "cnm = cnm.fit(images)" + "cnm.fit(images)" ] }, { @@ -466,7 +466,7 @@ " stride=stride, \n", " only_init_patch=True)\n", "cnm.params.set('spatial', {'se': np.ones((3,3,1), dtype=np.uint8)})\n", - "cnm = cnm.fit(images)\n", + "cnm.fit(images)\n", "print(('Number of components:' + str(cnm.estimates.A.shape[-1])))" ] }, diff --git a/demos/notebooks/demo_dendritic.ipynb b/demos/notebooks/demo_dendritic.ipynb index a8ae0191f..d7b00ffc2 100644 --- a/demos/notebooks/demo_dendritic.ipynb +++ b/demos/notebooks/demo_dendritic.ipynb @@ -329,7 +329,7 @@ "source": [ "%%time\n", "cnm = cnmf.CNMF(n_processes, params=opts, dview=dview)\n", - "cnm = cnm.fit(images)" + "cnm.fit(images)" ] }, { @@ -597,7 +597,7 @@ "outputs": [], "source": [ "cnm_patches = cnmf.CNMF(n_processes, params=opts, dview=dview)\n", - "cnm_patches = cnm_patches.fit(images)" + "cnm_patches.fit(images)" ] }, { diff --git a/demos/notebooks/demo_pipeline.ipynb b/demos/notebooks/demo_pipeline.ipynb index 42d58d441..e4baa3955 100644 --- a/demos/notebooks/demo_pipeline.ipynb +++ b/demos/notebooks/demo_pipeline.ipynb @@ -815,7 +815,7 @@ "outputs": [], "source": [ "%%time\n", - "cnmf_fit = cnmf_model.fit(images)" + "cnmf_model.fit(images)" ] }, { @@ -833,7 +833,7 @@ "metadata": {}, "outputs": [], "source": [ - "cnmf_fit.estimates.plot_contours_nb(img=correlation_image);" + "cnmf_model.estimates.plot_contours_nb(img=correlation_image);" ] }, { @@ -861,7 +861,7 @@ "outputs": [], "source": [ "%%time\n", - "cnmf_refit = cnmf_fit.refit(images, dview=cluster)" + "cnmf_refit = cnmf_model.refit(images, dview=cluster)" ] }, { diff --git a/demos/notebooks/demo_pipeline_voltage_imaging.ipynb b/demos/notebooks/demo_pipeline_voltage_imaging.ipynb index d2fcbea36..66095fd21 100644 --- a/demos/notebooks/demo_pipeline_voltage_imaging.ipynb +++ b/demos/notebooks/demo_pipeline_voltage_imaging.ipynb @@ -94,9 +94,9 @@ "outputs": [], "source": [ "# File path to movie file (will download if not present)\n", - "fnames = download_demo('demo_voltage_imaging.hdf5', 'volpy') \n", + "fnames = download_demo('demo_voltage_imaging.hdf5') \n", "# File path to ROIs file (will download if not present)\n", - "path_ROIs = download_demo('demo_voltage_imaging_ROIs.hdf5', 'volpy') \n", + "path_ROIs = download_demo('demo_voltage_imaging_ROIs.hdf5') \n", "file_dir = os.path.split(fnames)[0]\n", "\n", "# Setup some parameters for data and motion correction dataset parameters\n", diff --git a/demos/notebooks/demo_seeded_CNMF.ipynb b/demos/notebooks/demo_seeded_CNMF.ipynb index 17ec16671..fe3b3e8ae 100644 --- a/demos/notebooks/demo_seeded_CNMF.ipynb +++ b/demos/notebooks/demo_seeded_CNMF.ipynb @@ -29,6 +29,7 @@ "import logging\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "import os\n", "import cv2\n", "\n", "from glob import glob\n", @@ -221,7 +222,7 @@ "metadata": {}, "outputs": [], "source": [ - "crd = nb_plot_contour(mR, Ain.astype('float32'), mR.shape[0], mR.shape[1], thr=0.99)" + "nb_plot_contour(mR, Ain.astype('float32'), mR.shape[0], mR.shape[1], thr=0.99)" ] }, { @@ -445,7 +446,7 @@ "# For comparison, we can also run the unseeded CNMF. For this we have to change the parameters rf and only_init\n", "opts.change_params({'rf': 48, 'K': 12, 'merge_parallel': True})\n", "cnm = cnmf.CNMF(n_processes, params= opts, dview=dview)\n", - "cnm = cnm.fit(images)" + "cnm.fit(images)" ] }, { @@ -561,7 +562,7 @@ "Ain, mR = cm.base.rois.extract_binary_masks_from_structural_channel(\n", " cm.load(fname), expand_method='dilation', selem=np.ones((1, 1)))\n", "plt.figure()\n", - "crd = cm.utils.visualization.plot_contours(Ain.astype('float32'), mR, thr=0.99, display_numbers=False)\n", + "cm.utils.visualization.plot_contours(Ain.astype('float32'), mR, thr=0.99, display_numbers=False)\n", "plt.title('Contour plots of detected ROIs in the structural channel')\n", "'''" ] @@ -660,7 +661,8 @@ "plt.rcParams['figure.figsize'] = [16, 10]\n", "\n", "# load binary mask\n", - "mask = np.asarray(imread('/Users/hheiser/Desktop/testing data/file_00020_no_motion/avg_mask_fixed.png'), dtype=bool)\n", + "maskfile = os.path.join(cm.paths.caiman_datadir(), 'example_movies', 'avg_mask_fixed.png')\n", + "mask = np.asarray(imread(maskfile), dtype=bool)\n", "\n", "# calculate distances from nearest edge\n", "distances = ndi.distance_transform_edt(mask)\n", diff --git a/docs/source/core_functions.rst b/docs/source/core_functions.rst index 691f9c184..173f87fd1 100644 --- a/docs/source/core_functions.rst +++ b/docs/source/core_functions.rst @@ -121,8 +121,6 @@ CNMF .. automethod:: CNMF.update_temporal .. automethod:: CNMF.compute_residuals .. automethod:: CNMF.remove_components -.. automethod:: CNMF.HALS4traces -.. automethod:: CNMF.HALS4footprints .. automethod:: CNMF.merge_comps .. automethod:: CNMF.initialize .. automethod:: CNMF.preprocess diff --git a/environment-minimal.yml b/environment-minimal.yml index 81fa6898c..d7e9b3647 100644 --- a/environment-minimal.yml +++ b/environment-minimal.yml @@ -11,6 +11,7 @@ dependencies: - ipywidgets - matplotlib - moviepy +- pyside6 - pytest - numpy <2.0.0,>=1.26 - numpydoc diff --git a/environment.yml b/environment.yml index 014c3532c..eaa692baf 100644 --- a/environment.yml +++ b/environment.yml @@ -27,6 +27,7 @@ dependencies: - psutil - pynwb - pyqtgraph +- pyside6 - scikit-image >=0.19.0 - scikit-learn >=1.2 - scipy >= 1.10.1 diff --git a/example_movies/avg_mask_fixed.png b/example_movies/avg_mask_fixed.png new file mode 100644 index 000000000..4c0a12e86 Binary files /dev/null and b/example_movies/avg_mask_fixed.png differ diff --git a/pyproject.toml b/pyproject.toml index 96941901d..853b44aea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,34 +8,40 @@ authors = [ { name = "Valentina Staneva" }, { name = "Ben Deverett" }, { name = "Erick Cobos" }, - { name = "Jeremie Kalfon"}] + { name = "Jeremie Kalfon"}, + { name = "Kushal Kolar"}] readme = "README.md" dynamic = ["classifiers", "keywords", "license", "scripts", "version"] dependencies = [ "av", - "bokeh", + "bokeh >= 3.1.1", "coverage", "cython", - "h5py", - "holoviews", + "h5py >= 3.4.0", + "holoviews >= 1.16.2", "ipykernel", "ipython", "ipyparallel", "ipywidgets", "keras", "matplotlib", + "moviepy", "mypy", - "numpy", + "numpy <2.0.0, >=1.26", "numpydoc", "opencv-python", - "peakutils", + "panel >= 1.0.2", + "peakutils >= 1.3.5", "pims", "psutil", "pynwb", - "scikit-image", - "scikit-learn", - "scipy", - "tensorflow", + "pyside6", + "pytest", + "pytest-cov", + "scikit-image >= 0.19.0", + "scikit-learn >= 1.2", + "scipy >= 1.10.1", + "tensorflow >= 2.4.0, <2.16", "tifffile", "tqdm", "yapf", @@ -45,6 +51,7 @@ dependencies = [ [project.optional-dependencies] jupyter = [ "jupyter", + "jupyter_bokeh", "pyqtgraph", "tk" ] diff --git a/setup.py b/setup.py index 2cb65e779..02845a582 100755 --- a/setup.py +++ b/setup.py @@ -33,8 +33,8 @@ extra_dirs = ['bin', 'demos', 'docs', 'model'] data_files = [('share/caiman', ['LICENSE.txt', 'README.md', 'test_demos.sh', 'VERSION']), - ('share/caiman/example_movies', ['example_movies/data_endoscope.tif', 'example_movies/demoMovie.tif']), - ('share/caiman/testdata', ['testdata/groundtruth.npz', 'testdata/example.npz', 'testdata/2d_sbx.mat', 'testdata/2d_sbx.sbx', 'testdata/3d_sbx_1.mat', 'testdata/3d_sbx_1.sbx', 'testdata/3d_sbx_2.mat', 'testdata/3d_sbx_2.sbx']), + ('share/caiman/example_movies', ['example_movies/data_endoscope.tif', 'example_movies/demoMovie.tif', 'example_movies/avg_mask_fixed.png']), + ('share/caiman/testdata', ['testdata/groundtruth.npz', 'testdata/example.npz', 'testdata/2d_sbx.mat', 'testdata/2d_sbx.sbx', 'testdata/2d_sbx_bidi.mat', 'testdata/2d_sbx_bidi.sbx', 'testdata/3d_sbx_1.mat', 'testdata/3d_sbx_1.sbx', 'testdata/3d_sbx_2.mat', 'testdata/3d_sbx_2.sbx']), ] for part in extra_dirs: newpart = [("share/caiman/" + d, [os.path.join(d,f) for f in files]) for d, folders, files in os.walk(part)] diff --git a/testdata/2d_sbx_bidi.mat b/testdata/2d_sbx_bidi.mat new file mode 100644 index 000000000..9eb119cda Binary files /dev/null and b/testdata/2d_sbx_bidi.mat differ diff --git a/testdata/2d_sbx_bidi.sbx b/testdata/2d_sbx_bidi.sbx new file mode 100644 index 000000000..dfb220fe6 Binary files /dev/null and b/testdata/2d_sbx_bidi.sbx differ