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..8211f5549 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. @@ -163,7 +168,6 @@ Special thanks to the following people for letting us use their datasets in demo Also a special thanks to: * Eric Thompson, for various strong contributions to code and demos, both before and during his employment at the Flatiron Institute. -* Cai Changjia, for Volpy * Ethan Blackwood, for several contributions in various areas # License diff --git a/VERSION b/VERSION index e6dbb7c23..6b89d58f8 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.11.5 +1.12.2 diff --git a/bin/caiman_gui.py b/bin/caiman_gui.py index 0a02eb459..cd54a4216 100755 --- a/bin/caiman_gui.py +++ b/bin/caiman_gui.py @@ -16,6 +16,23 @@ from caiman.source_extraction.cnmf.cnmf import load_CNMF from caiman.source_extraction.cnmf.params import CNMFParams +############################## +# This is a tool that can help you visualise dumps of the CNMF object, +# generated as demonstrated in demo_pipeline. +# +# If you're already using it, that's fine, but there are better tools +# out there (such as the Mesmerize-vis package). To use, point it at the +# hdf5 file that CLI demo or jupyter notebook made, or have your own code +# call the save method on the CNMF object, and then feed that (and then the +# mmap file) into this script when it asks, and you'll see your ROIs and +# be able to experiment with thresholding (and re-save with those different +# thresholds). +# +# mesmerize-vis offers a better approach to these things, and is much more +# recent and better maintained. + + + # Always start by initializing Qt (only once per application) app = QtWidgets.QApplication([]) @@ -66,7 +83,7 @@ def make_color_img(img, gain=255, min_max=None, out_type=np.uint8): M = FileDialog() d, f = os.path.split(cnm_obj.mmap_file) cnm_obj.mmap_file = M.getOpenFileName(caption='Load memory mapped file', - directory=d, filter=f + ';;*.mmap')[0] + filter=f + ';;*.mmap')[0] if fpath[-3:] == 'nwb': mov = caiman.load(cnm_obj.mmap_file, @@ -405,15 +422,15 @@ def move(event): ## PARAMS -params = [{'name': 'min_cnn_thr', 'type': 'float', 'value': 0.99, 'limits': (0, 1),'step':0.01}, - {'name': 'cnn_lowest', 'type': 'float', 'value': 0.1, 'limits': (0, 1),'step':0.01}, - {'name': 'rval_thr', 'type': 'float', 'value': 0.85, 'limits': (-1, 1),'step':0.01}, - {'name': 'rval_lowest', 'type': 'float', 'value': -1, 'limits': (-1, 1),'step':0.01}, - {'name': 'min_SNR', 'type': 'float', 'value': 2, 'limits': (0, 20),'step':0.1}, - {'name': 'SNR_lowest', 'type': 'float', 'value': 0, 'limits': (0, 20),'step':0.1}, - {'name': 'RESET', 'type': 'action'}, +params = [{'name': 'min_cnn_thr', 'type': 'float', 'value': 0.99, 'limits': (0, 1), 'step':0.01}, + {'name': 'cnn_lowest', 'type': 'float', 'value': 0.1, 'limits': (0, 1), 'step':0.01}, + {'name': 'rval_thr', 'type': 'float', 'value': 0.85, 'limits': (-1, 1), 'step':0.01}, + {'name': 'rval_lowest', 'type': 'float', 'value': -1, 'limits': (-1, 1), 'step':0.01}, + {'name': 'min_SNR', 'type': 'float', 'value': 2, 'limits': (0, 20), 'step':0.1}, + {'name': 'SNR_lowest', 'type': 'float', 'value': 0, 'limits': (0, 20), 'step':0.1}, + {'name': 'RESET', 'type': 'action'}, {'name': 'SHOW BACKGROUND', 'type': 'action'}, - {'name': 'SHOW NEURONS', 'type': 'action'} + {'name': 'SHOW NEURONS', 'type': 'action'} ] ## Create tree of Parameter objects @@ -421,7 +438,7 @@ def move(event): params_action = [{'name': 'Filter components', 'type': 'bool', 'value': True, 'tip': "Filter components"}, - {'name': 'View components', 'type': 'list', 'values': ['All','Accepted', + {'name': 'View components', 'type': 'list', 'limits': ['All','Accepted', 'Rejected', 'Unassigned'], 'value': 'All'}, {'name': 'ADD GROUP', 'type': 'action'}, {'name': 'REMOVE GROUP', 'type': 'action'}, diff --git a/caiman/__init__.py b/caiman/__init__.py index 7a5a04343..a2b3416b1 100644 --- a/caiman/__init__.py +++ b/caiman/__init__.py @@ -1,10 +1,11 @@ #!/usr/bin/env python -import pkg_resources +import importlib.metadata + from caiman.base.movies import movie, load, load_movie_chain, _load_behavior, play_movie from caiman.base.timeseries import concatenate from caiman.cluster import start_server, stop_server from caiman.mmapping import load_memmap, save_memmap, save_memmap_each, save_memmap_join from caiman.summary_images import local_correlations -__version__ = pkg_resources.get_distribution('caiman').version +__version__ = importlib.metadata.version('caiman') diff --git a/caiman/base/movies.py b/caiman/base/movies.py index f0a9a5796..40676ec05 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 @@ -271,7 +187,7 @@ def bin_median(self, window: int = 10) -> np.ndarray: median image """ - T, d1, d2 = np.shape(self) + T, d1, d2 = self.shape num_windows = int(T // window) num_frames = num_windows * window return np.nanmedian(np.nanmean(np.reshape(self[:num_frames], (window, num_windows, d1, d2)), axis=0), axis=0) @@ -291,7 +207,7 @@ def bin_median_3d(self, window=10): median image """ - T, d1, d2, d3 = np.shape(self) + T, d1, d2, d3 = self.shape num_windows = int(T // window) num_frames = num_windows * window return np.nanmedian(np.nanmean(np.reshape(self[:num_frames], (window, num_windows, d1, d2, d3)), axis=0), @@ -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 @@ -598,18 +465,18 @@ def computeDFF(self, secsWindow: int = 5, quantilMin: int = 8, method: str = 'on if np.min(self) <= 0 and method != 'only_baseline': raise ValueError("All pixels must be positive") - numFrames, linePerFrame, pixPerLine = np.shape(self) + numFrames, linePerFrame, pixPerLine = self.shape downsampfact = int(secsWindow * self.fr) logger.debug(f"Downsample factor: {downsampfact}") elm_missing = int(np.ceil(numFrames * 1.0 / downsampfact) * downsampfact - numFrames) padbefore = int(np.floor(elm_missing / 2.)) padafter = int(np.ceil(elm_missing / 2.)) - logger.debug(f'Initial Size Image: {np.shape(self)}') + logger.debug(f'Initial Size Image: {self.shape}') sys.stdout.flush() mov_out = movie(np.pad(self.astype(np.float32), ((padbefore, padafter), (0, 0), (0, 0)), mode='reflect'), **self.__dict__) - numFramesNew, linePerFrame, pixPerLine = np.shape(mov_out) + numFramesNew, linePerFrame, pixPerLine = mov_out.shape # compute baseline quickly logger.debug("binning data ...") @@ -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 @@ -697,7 +540,7 @@ def IPCA(self, components: int = 50, batch: int = 1000) -> tuple[np.ndarray, np. """ # vectorize the images - num_frames, h, w = np.shape(self) + num_frames, h, w = self.shape frame_size = h * w frame_samples = np.reshape(self, (num_frames, frame_size)).T @@ -768,7 +611,7 @@ def IPCA_stICA(self, joint_ics = ica.fit_transform(eigenstuff) # extract the independent frames - _, h, w = np.shape(self) + _, h, w = self.shape frame_size = h * w ind_frames = joint_ics[:frame_size, :] ind_frames = np.reshape(ind_frames.T, (componentsICA, h, w)) @@ -780,7 +623,7 @@ def IPCA_denoise(self, components: int = 50, batch: int = 1000): Create a denoised version of the movie using only the first 'components' components """ _, _, clean_vectors = self.IPCA(components, batch) - self = self.__class__(np.reshape(np.float32(clean_vectors.T), np.shape(self)), **self.__dict__) + self = self.__class__(np.reshape(np.float32(clean_vectors.T), self.shape), **self.__dict__) return self def local_correlations(self, @@ -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, @@ -1451,7 +1214,7 @@ def load(file_name: Union[str, list[str]], if input_arr.ndim == 2: if shape is not None: - _, T = np.shape(input_arr) + _, T = input_arr.shape d1, d2 = shape input_arr = np.transpose(np.reshape(input_arr, (d1, d2, T), order='F'), (2, 0, 1)) else: @@ -1598,13 +1361,13 @@ def load_movie_chain(file_list: list[str], if m.ndim == 2: m = m[np.newaxis, :, :] - _, h, w = np.shape(m) + _, h, w = m.shape m = m[:, top:h - bottom, left:w - right] else: if m.ndim == 3: m = m[np.newaxis, :, :, :] - _, h, w, d = np.shape(m) + _, h, w, d = m.shape m = m[:, top:h - bottom, left:w - right, z_top:d - z_bottom] mov.append(m) @@ -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/base/rois.py b/caiman/base/rois.py index d1c11c374..186918516 100644 --- a/caiman/base/rois.py +++ b/caiman/base/rois.py @@ -123,14 +123,14 @@ def extract_binary_masks_from_structural_channel(Y, def mask_to_2d(mask): # todo todocument if mask.ndim > 2: - _, d1, d2 = np.shape(mask) + _, d1, d2 = mask.shape dims = d1, d2 return scipy.sparse.coo_matrix(np.reshape(mask[:].transpose([1, 2, 0]), ( np.prod(dims), -1, ), order='F')) else: - dims = np.shape(mask) + dims = mask.shape return scipy.sparse.coo_matrix(np.reshape(mask, ( np.prod(dims), -1, @@ -140,7 +140,7 @@ def mask_to_2d(mask): def get_distance_from_A(masks_gt, masks_comp, min_dist=10) -> list: # todo todocument - _, d1, d2 = np.shape(masks_gt) + _, d1, d2 = masks_gt.shape dims = d1, d2 A_ben = scipy.sparse.csc_matrix(np.reshape(masks_gt[:].transpose([1, 2, 0]), ( np.prod(dims), @@ -214,7 +214,7 @@ def nf_match_neurons_in_binary_masks(masks_gt, """ logger = logging.getLogger("caiman") - _, d1, d2 = np.shape(masks_gt) + _, d1, d2 = masks_gt.shape dims = d1, d2 # transpose to have a sparse list of components, then reshaping it to have a 1D matrix red in the Fortran style @@ -244,8 +244,8 @@ def nf_match_neurons_in_binary_masks(masks_gt, # compute precision and recall TP = np.sum(np.array(costs) < thresh_cost) * 1. - FN = np.shape(masks_gt)[0] - TP - FP = np.shape(masks_comp)[0] - TP + FN = masks_gt.shape[0] - TP + FP = masks_comp.shape[0] - TP TN = 0 performance = dict() @@ -259,9 +259,9 @@ def nf_match_neurons_in_binary_masks(masks_gt, idx_tp_ben = matches[0][idx_tp] # ground truth idx_tp_cnmf = matches[1][idx_tp] # algorithm - comp - idx_fn = np.setdiff1d(list(range(np.shape(masks_gt)[0])), matches[0][idx_tp]) + idx_fn = np.setdiff1d(list(range(masks_gt.shape[0])), matches[0][idx_tp]) - idx_fp = np.setdiff1d(list(range(np.shape(masks_comp)[0])), matches[1][idx_tp]) + idx_fp = np.setdiff1d(list(range(masks_comp.shape[0])), matches[1][idx_tp]) idx_fp_cnmf = idx_fp @@ -736,8 +736,8 @@ def distance_masks(M_s:list, cm_s: list[list], max_dist: float, enclosed_thr: Op test_comp = test_comp.copy()[:, :] # the number of components for each - nb_gt = np.shape(gt_comp)[-1] - nb_test = np.shape(test_comp)[-1] + nb_gt = gt_comp.shape[-1] + nb_test = test_comp.shape[-1] D = np.ones((nb_gt, nb_test)) cmgt_comp = np.array(cmgt_comp) diff --git a/caiman/behavior/behavior.py b/caiman/behavior/behavior.py index 9863f5603..decbcbef7 100644 --- a/caiman/behavior/behavior.py +++ b/caiman/behavior/behavior.py @@ -46,7 +46,7 @@ def select_roi(img: np.ndarray, n_rois: int = 1) -> list: fig = plt.figure() plt.imshow(img, cmap=matplotlib.cm.gray) pts = fig.ginput(0, timeout=0) - mask = np.zeros(np.shape(img), dtype=np.int32) + mask = np.zeros(img.shape, dtype=np.int32) pts = np.asarray(pts, dtype=np.int32) cv2.fillConvexPoly(mask, pts, (1, 1, 1), lineType=cv2.LINE_AA) masks.append(mask) @@ -317,11 +317,11 @@ def extract_components(mov_tot, mov_tot = mov_tot / norm_fact[:, np.newaxis, np.newaxis, np.newaxis] else: norm_fact = np.array([1., 1.]) - c, T, d1, d2 = np.shape(mov_tot) + c, T, d1, d2 = mov_tot.shape else: norm_fact = 1 - T, d1, d2 = np.shape(mov_tot) + T, d1, d2 = mov_tot.shape c = 1 tt = time.time() diff --git a/caiman/caimanmanager.py b/caiman/caimanmanager.py index 968bbfca6..82a79e8ce 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') ] ############### @@ -126,7 +127,7 @@ def do_check_install(targdir: str, inplace: bool = False) -> None: def do_run_pytest(targdir: str) -> None: - out, err, ret = runcmd(["pytest", "--verbose", "--pyargs", "caiman"]) + out, err, ret = runcmd(["pytest", "--verbose", "-W", "ignore::sklearn.exceptions.ConvergenceWarning", "--pyargs", "caiman"]) if ret != 0: print(f"pytest failed with return code {ret}") sys.exit(ret) diff --git a/caiman/cluster.py b/caiman/cluster.py index dde1b336f..19bae96b6 100644 --- a/caiman/cluster.py +++ b/caiman/cluster.py @@ -222,12 +222,12 @@ def setup_cluster(backend:str = 'multiprocessing', if backend == 'multiprocessing' or backend == 'local': if backend == 'local': - logger.warn('The local backend is an alias for the multiprocessing backend, and the alias may be removed in some future version of Caiman') + logger.warning('The local backend is an alias for the multiprocessing backend, and the alias may be removed in some future version of Caiman') if len(multiprocessing.active_children()) > 0: if ignore_preexisting: - logger.warn('Found an existing multiprocessing pool. ' - 'This is often indicative of an already-running CaImAn cluster. ' - 'You have configured the cluster setup to not raise an exception.') + logger.warning('Found an existing multiprocessing pool. ' + 'This is often indicative of an already-running CaImAn cluster. ' + 'You have configured the cluster setup to not raise an exception.') else: raise Exception( 'A cluster is already running. Terminate with dview.terminate() if you want to restart.') @@ -256,7 +256,7 @@ def setup_cluster(backend:str = 'multiprocessing', elif backend == "single" or single_thread: if single_thread: - logger.warn('The single_thread flag to setup_cluster() is deprecated and may be removed in the future') + logger.warning('The single_thread flag to setup_cluster() is deprecated and may be removed in the future') dview = None c = None n_processes = 1 diff --git a/caiman/components_evaluation.py b/caiman/components_evaluation.py index 5bf9a47c0..e1622f3cc 100644 --- a/caiman/components_evaluation.py +++ b/caiman/components_evaluation.py @@ -67,14 +67,17 @@ def compute_event_exceptionality(traces: np.ndarray, erfc: ndarray probability at each time step of observing the N consecutive actual trace values given the distribution of noise - noise_est: ndarray - the components ordered according to the fitness + std_r: + Standard deviation of r + + mode: + Mode of the traces """ if N == 0: # Without this, numpy ranged syntax does not work correctly, and also N=0 is conceptually incoherent raise Exception("FATAL: N=0 is not a valid value for compute_event_exceptionality()") - T = np.shape(traces)[-1] + T = traces.shape[-1] if use_mode_fast: md = caiman.utils.stats.mode_robust_fast(traces, axis=1) else: @@ -85,7 +88,6 @@ def compute_event_exceptionality(traces: np.ndarray, # only consider values under the mode to determine the noise standard deviation ff1 = -ff1 * (ff1 < 0) if robust_std: - # compute 25 percentile ff1 = np.sort(ff1, axis=1) ff1[ff1 == 0] = np.nan @@ -138,7 +140,7 @@ def find_activity_intervals(C, Npeaks: int = 5, tB=-3, tA=10, thres: float = 0.3 # todo todocument logger = logging.getLogger("caiman") - K, T = np.shape(C) + K, T = C.shape L:list = [] for i in range(K): if np.sum(np.abs(np.diff(C[i, :]))) == 0: @@ -208,7 +210,7 @@ def classify_components_ep(Y, A, C, b, f, Athresh=0.1, Npeaks=5, tB=-3, tA=10, t """ logger = logging.getLogger("caiman") - K, _ = np.shape(C) + K, _ = C.shape A = csc_matrix(A) AA = (A.T * A).toarray() nA = np.sqrt(np.array(A.power(2).sum(0))) @@ -415,7 +417,7 @@ def evaluate_components(Y: np.ndarray, tB = np.minimum(-2, np.floor(-5. / 30 * final_frate)) tA = np.maximum(5, np.ceil(25. / 30 * final_frate)) logger.info(f'{tB=},{tA=}') - dims, T = np.shape(Y)[:-1], np.shape(Y)[-1] + dims, T = Y.shape[:-1], Y.shape[-1] Yr = np.reshape(Y, (np.prod(dims), T), order='F') @@ -427,7 +429,7 @@ def evaluate_components(Y: np.ndarray, logger.debug('Removing Baseline') if remove_baseline: - num_samps_bl = np.minimum(np.shape(traces)[-1]// 5, 800) + num_samps_bl = np.minimum(traces.shape[-1]// 5, 800) slow_baseline = False if slow_baseline: @@ -441,7 +443,7 @@ def evaluate_components(Y: np.ndarray, padbefore = int(np.floor(elm_missing / 2.)) padafter = int(np.ceil(elm_missing / 2.)) tr_tmp = np.pad(traces.T, ((padbefore, padafter), (0, 0)), mode='reflect') - numFramesNew, num_traces = np.shape(tr_tmp) + numFramesNew, num_traces = tr_tmp.shape # compute baseline quickly logger.debug("binning data ...") tr_BL = np.reshape(tr_tmp, (downsampfact, numFramesNew // downsampfact, num_traces), order='F') diff --git a/caiman/mmapping.py b/caiman/mmapping.py index a69558169..e52988fe0 100644 --- a/caiman/mmapping.py +++ b/caiman/mmapping.py @@ -577,7 +577,7 @@ def parallel_dot_product(A: np.ndarray, b, block_size: int = 5000, dview=None, t logger = logging.getLogger("caiman") pars = [] - d1, d2 = np.shape(A) + d1, d2 = A.shape b = pickle.dumps(b) logger.debug(f'parallel dot product block size: {block_size}') @@ -598,9 +598,9 @@ def parallel_dot_product(A: np.ndarray, b, block_size: int = 5000, dview=None, t b = pickle.loads(b) if transpose: - output = np.zeros((d2, np.shape(b)[-1]), dtype=np.float32) + output = np.zeros((d2, b.shape[-1]), dtype=np.float32) else: - output = np.zeros((d1, np.shape(b)[-1]), dtype=np.float32) + output = np.zeros((d1, b.shape[-1]), dtype=np.float32) if dview is None: if transpose: diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py index f5c030dbf..722407ae4 100644 --- a/caiman/motion_correction.py +++ b/caiman/motion_correction.py @@ -202,7 +202,22 @@ def __init__(self, fname, min_mov=None, dview=None, max_shifts=(6, 6), niter_rig self.indices = indices self.shifts_interpolate = shifts_interpolate if self.use_cuda: - logger.warn("cuda is no longer supported; this kwarg will be removed in a future version of caiman") + logger.warning("cuda is no longer supported; this kwarg will be removed in a future version of caiman") + + def __str__(self): + ret = f"Caiman MotionCorrect Object. Subfields:{list(self.__dict__.keys())}" + if hasattr(self, 'fname'): + ret += "MotionCorrect backing fname is {self.fname}" + if 'tmp_mov_mot_corr' in self.fname: + ret += "Target data was passed inline" + return ret + + def __repr__(self): + return self.__str__() + + def __getitem__(self, idx): + return getattr(self, idx) + # We intentionally do not provide __setitem__ def motion_correct(self, template=None, save_movie=False): """general function for performing all types of motion correction. The @@ -218,12 +233,8 @@ def motion_correct(self, template=None, save_movie=False): save_movie: bool, default: False flag for saving motion corrected file(s) as memory mapped file(s) - - Returns: - self """ - # TODO: Review the docs here, and also why we would ever return self - # from a method that is not a constructor + # TODO: Review the docs here if self.min_mov is None: if self.gSig_filt is None: iterator = caiman.base.movies.load_iter(self.fname[0], @@ -255,7 +266,6 @@ def motion_correct(self, template=None, save_movie=False): b0 = np.ceil(np.max(np.abs(self.shifts_rig))) self.border_to_0 = b0.astype(int) self.mmap_file = self.fname_tot_els if self.pw_rigid else self.fname_tot_rig - return self def motion_correct_rigid(self, template: Optional[np.ndarray] = None, save_movie=False) -> None: """ @@ -483,7 +493,8 @@ def apply_shifts_movie(self, fname, rigid_shifts: Optional[bool] = None, save_me return caiman.movie(m_reg) def apply_shift_iteration(img, shift, border_nan=False, border_type=cv2.BORDER_REFLECT): - # todo todocument + # Used by MotionCorrect.apply_shifts_movie(), tile_and_correct(), and apply_shift_online() + # Applies an xy shift to one iterable unit of a movie. sh_x_n, sh_y_n = shift w_i, h_i = img.shape @@ -717,7 +728,7 @@ def motion_correct_online_multifile(list_files, add_to_movie, order='C', **kwarg kwargs_['order'] = order total_frames = 0 for file_ in list_files: - logger.info('Processing: {file_}') + logger.info(f'Processing: {file_}') kwargs_['template'] = template kwargs_['save_base_name'] = os.path.splitext(file_)[0] tffl = tifffile.TiffFile(file_) @@ -992,7 +1003,7 @@ def bin_median(mat, window=10, exclude_nans=True): median image """ - T, d1, d2 = np.shape(mat) + T, d1, d2 = mat.shape if T < window: window = T num_windows = int(T // window) @@ -1021,7 +1032,7 @@ def bin_median_3d(mat, window=10, exclude_nans=True): median image """ - T, d1, d2, d3 = np.shape(mat) + T, d1, d2, d3 = mat.shape if T < window: window = T num_windows = int(T // window) @@ -1281,6 +1292,32 @@ def _compute_error(cross_correlation_max, src_amp, target_amp): (src_amp * target_amp) return np.sqrt(np.abs(error)) +def init_cuda_process(): + """ + Initialize a PyCUDA context at global scope so that it can be accessed + from processes when using multithreading + """ + global cudactx + + cudadrv.init() + dev = cudadrv.Device(0) + cudactx = dev.make_context() # type: ignore + atexit.register(cudactx.pop) # type: ignore + + +def close_cuda_process(n): + """ + Cleanup cuda process + """ + + global cudactx + + import skcuda.misc as cudamisc + try: + cudamisc.done_context(cudactx) # type: ignore + except: + pass + def register_translation_3d(src_image, target_image, upsample_factor = 1, space = "real", shifts_lb = None, shifts_ub = None, max_shifts = [10,10,10]): @@ -1682,14 +1719,14 @@ def apply_shifts_dft(src_freq, shifts, diffphase, is_freq=True, border_nan=True) src_freq = np.array(src_freq, dtype=np.complex128, copy=False) if not is3D: - nr, nc = np.shape(src_freq) + nr, nc = src_freq.shape Nr = ifftshift(np.arange(-np.fix(nr/2.), np.ceil(nr/2.))) Nc = ifftshift(np.arange(-np.fix(nc/2.), np.ceil(nc/2.))) Nc, Nr = np.meshgrid(Nc, Nr) Greg = src_freq * np.exp(1j * 2 * np.pi * (-shifts[0] * Nr / nr - shifts[1] * Nc / nc)) else: - nr, nc, nd = np.array(np.shape(src_freq), dtype=float) + nr, nc, nd = np.array(src_freq.shape, dtype=float) Nr = ifftshift(np.arange(-np.fix(nr / 2.), np.ceil(nr / 2.))) Nc = ifftshift(np.arange(-np.fix(nc / 2.), np.ceil(nc / 2.))) Nd = ifftshift(np.arange(-np.fix(nd / 2.), np.ceil(nd / 2.))) @@ -2639,10 +2676,10 @@ def compute_metrics_motion_correction(fname, final_size_x, final_size_y, swap_di m_min = mi + (ma - mi) / 100 m_max = mi + (ma - mi) / 4 - max_shft_x = int(np.ceil((np.shape(m)[1] - final_size_x) / 2)) - max_shft_y = int(np.ceil((np.shape(m)[2] - final_size_y) / 2)) - max_shft_x_1 = - ((np.shape(m)[1] - max_shft_x) - (final_size_x)) - max_shft_y_1 = - ((np.shape(m)[2] - max_shft_y) - (final_size_y)) + max_shft_x = int(np.ceil((m.shape[1] - final_size_x) / 2)) + max_shft_y = int(np.ceil((m.shape[2] - final_size_y) / 2)) + max_shft_x_1 = - ((m.shape[1] - max_shft_x) - (final_size_x)) + max_shft_y_1 = - ((m.shape[2] - max_shft_y) - (final_size_y)) if max_shft_x_1 == 0: max_shft_x_1 = None diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 5143ab6db..a55ca93b1 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -24,6 +24,7 @@ import shutil import pathlib import psutil +import pynwb import scipy import sys @@ -89,133 +90,174 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p sniper_mode=False, use_peak_max=False, test_both=False, expected_comps=500, params=None): """ - Constructor of the CNMF method + Constructor of CNMF objects + + Please pass most or all params in as a Params object, using the params keyword. Most of the other + arguments are legacy and may be removed in some (likely distant) future update. + + If you use the params argument, most of the other arguments will be ignored. + These arguments are exceptions: dview, skip_refinement, remove_very_bad_comps, Ain, Cin, b_in, f_in Args: n_processes: int - number of processed used (if in parallel this controls memory usage) + number of processes used (if in parallel this controls memory usage) k: int - number of neurons expected per FOV (or per patch if patches_pars is None) + number of neurons expected per FOV (or per patch if patches_pars is None) gSig: tuple expected half size of neurons + gSiz: tuple + half-size of bounding box for each neuron (computed from gSig is None) + merge_thresh: float merging threshold, max correlation allowed - dview: Direct View object - for parallelization purposes when using ipyparallel - p: int order of the autoregressive process used to estimate deconvolution + dview: Direct View object + for parallelization purposes when using ipyparallel + Ain: np.ndarray - if known, it is the initial estimate of spatial filters. Array must be of type `bool` in 'F' order of shape: [n_pixels, n_components] + if known, it is the initial estimate of spatial filters. Array must be of type `bool` in 'F' order of shape: [n_pixels, n_components]. + Used to build estimates + + Cin - Used to build estimates + + b_in - Used to build estimates + + f_in - Used to build estimates + + do_merge: boolean + Whether to merge components in refinement ssub: int - downsampleing factor in space + downsampling factor in space tsub: int - downsampling factor in time + downsampling factor in time p_ssub: int downsampling factor in space for patches + p_tsub: int + downsampling factor in time for patches + method_init: str can be greedy_roi or sparse_nmf alpha_snmf: float weight of the sparsity regularization - p_tsub: int - downsampling factor in time for patches - - rf: int + rf:int half-size of the patches in pixels. rf=25, patches are 50x50 - gnb: int - number of global background components - - nb_patch: int - number of background components per patch - stride: int amount of overlap between the patches in pixels memory_fact: float unitless number accounting how much memory should be used. You will - need to try different values to see which one would work the default is OK for a 16 GB system + need to try different values to find a workable value. The default is OK for a 16 GB system + + gnb: int + number of global background components - N_samples_fitness: int - number of samples over which exceptional events are computed (See utilities.evaluate_components) + nb_patch: int + number of background components per patch - only_init_patch= boolean + only_init_patch: boolean only run initialization on patches + only_init_patch (not documented) + method_deconvolution = 'oasis' or 'cvxpy' method used for deconvolution. Suggested 'oasis' see Friedrich J, Zhou P, Paninski L. Fast Online Deconvolution of Calcium Imaging Data. PLoS Comput Biol. 2017; 13(3):e1005423. - n_pixels_per_process: int. + n_pixels_per_process: int Number of pixels to be processed in parallel per core (no patch mode). Decrease if memory problems - block_size: int. + block_size_temp: int Number of pixels to be used to perform residual computation in blocks. Decrease if memory problems num_blocks_per_run_spat: int - In case of memory problems you can reduce this numbers, controlling the number of blocks processed in parallel during residual computing + In case of memory problems you can reduce this number, reducing the number of blocks processed in parallel during residual computing num_blocks_per_run_temp: int - In case of memory problems you can reduce this numbers, controlling the number of blocks processed in parallel during residual computing + In case of memory problems you can reduce this number, reducing the number of blocks processed in parallel during residual computing - check_nan: Boolean. + check_nan: boolean Check if file contains NaNs (costly for very large files so could be turned off) - skip_refinement: - Bool. If true it only performs one iteration of update spatial update temporal instead of two + skip_refinement: boolean + If true it only performs one iteration of update spatial update temporal instead of two - normalize_init=Bool. - Differences in intensities on the FOV might caus troubles in the initialization when patches are not used, - so each pixels can be normalized by its median intensity + normalize_init: boolean + Differences in intensities in the FOV can cause troubles in the initialization when patches are not used. + This compensates by normalizing pixels by median intensity - options_local_NMF: + options_local_NMF: dict experimental, not to be used - remove_very_bad_comps:Bool - whether to remove components with very low values of component quality directly on the patch. - This might create some minor imprecisions. + minibatch_shape (not documented) + + minibatch_suff_strat (not documented) + + update_num_comps (not documented) + + rval_thr (not documented) + + thresh_fitness_delta (not documented) + + thresh_fitness_raw (not documented) + + thresh_overlap (not documented) + + max_comp_update_shape: + threshold number of components after which selective updating starts (using the parameter num_times_comp_updated) + + num_times_comp_updated: + number of times each component is updated. In inf components are updated at every initbatch time steps + + batch_update_suff_stat (not documented) + + s_min (not documented) + + remove_very_bad_comps: boolean + Whether to remove components with very low values of component quality directly on the patch. + This might create some minor imprecisions. However benefits can be considerable if done because if many components (>2000) are created and joined together, operation that causes a bottleneck - border_pix:int + border_pix: int number of pixels to not consider in the borders - low_rank_background:bool - if True the background is approximated with gnb components. If false every patch keeps its background (overlaps are randomly assigned to one spatial component only) - In the False case all the nonzero elements of the background components are updated using hals (to be used with one background per patch) + low_rank_background: boolean + Ff True the background is approximated with gnb components. If false every patch keeps its background (overlaps are randomly assigned to one spatial component only) + In the False case all the nonzero elements of the background components are updated using hals (to be used with one background per patch) - update_background_components:bool + update_background_components: boolean whether to update the background components during the spatial phase + rolling_sum (not documented) + + rolling_length (not documented) + min_corr: float minimal correlation peak for 1-photon imaging initialization min_pnr: float - minimal peak to noise ratio for 1-photon imaging initialization + minimal peak to noise ratio for 1-photon imaging initialization ring_size_factor: float - it's the ratio between the ring radius and neuron diameters. + Ratio between the ring radius and neuron diameters - max_comp_update_shape: - threshold number of components after which selective updating starts (using the parameter num_times_comp_updated) - - num_times_comp_updated: - number of times each component is updated. In inf components are updated at every initbatch time steps + center_psf (not documented) - expected_comps: int - number of expected components (try to exceed the expected) + use_dense (not documented) deconv_flag : bool, optional If True, deconvolution is also performed using OASIS @@ -225,14 +267,14 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p n_refit : int, optional Number of pools (cf. oasis.pyx) prior to the last one that are refitted when - simultaneously demixing and denoising/deconvolving. - - N_samples_exceptionality : int, optional - Number of consecutives intervals to be considered when testing new neuron candidates + simultaneously demixing and denoising/deconvolving del_duplicates: bool whether to delete the duplicated created in initialization + N_samples_exceptionality : int, optional + Number of consecutives intervals to be considered when testing new neuron candidates + max_num_added : int, optional maximum number of components to be added at each step in OnACID @@ -242,12 +284,31 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p thresh_CNN_noisy: float threshold on the per patch CNN classifier for online algorithm + fr: int + frame rate + + decay_time (not documented) + + min_SNR (not documented) + ssub_B: int, optional - downsampleing factor for 1-photon imaging background computation + downsampling factor for 1-photon imaging background computation init_iter: int, optional number of iterations for 1-photon imaging initialization + sniper_mode: boolean (undocumented) + + use_peak_max: boolean (undocumented) + + test_both: boolean (undocumented) + + expected_comps: int + number of expected components (try to exceed the expected) + + params: optional, CNMFParams object + If you have an object with all the params you want ready, you can pass it in rather than initialise in pieces as above + """ self.dview = dview @@ -291,9 +352,37 @@ 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 __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 + Packages the analysis pipeline (motion correction, memory mapping, patch based CNMF processing and component evaluation) in a single method that can be called on a specific (sequence of) file(s). It is assumed that the CNMF object already contains a params object @@ -378,12 +467,14 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False, outpu Cn[np.isnan(Cn)] = 0 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 + #fit_cnm.params.change_params({'p': self.params.get('preprocess', 'p')}) # RE-RUN 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 @@ -396,12 +487,24 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False, outpu for log_file in log_files: os.remove(log_file) + + if output_dir is not None: + # copy the result files to the specified output directory + output_dir = pathlib.Path(output_dir) + files_to_move = [fname_new, fname_init_hdf5, fname_hdf5] + if motion_correct: + files_to_move += fname_mc + for f in files_to_move: + f = pathlib.Path(f) + shutil.copy2(f, output_dir) + # 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 @@ -409,7 +512,7 @@ def fit_file(self, motion_correct=False, indices=None, include_eval=False, outpu def refit(self, images, dview=None): """ - Refits the data using CNMF initialized from a previous iteration + Refit data using CNMF initialized from a previous iteration Args: images @@ -426,20 +529,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 """ @@ -447,23 +549,25 @@ 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) if self.params.get('patch', 'rf') is None and (is_sliced or 'ndarray' in str(type(images))): images = images[tuple(indices)] self.dview = None - logger.info("Parallel processing in a single patch is not available for loaded in memory or sliced data.") + logger.info("Parallel processing in a single patch is not available for data that is in memory or sliced") 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): @@ -471,26 +575,20 @@ 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 + # FIXME The code below is really ugly and it's hard to tell if it's doing the right thing. + # These decisions should also probably be set higher up the call stack in some kind of a performance + # API (if we go with execution contexts, definitely there) 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'))) @@ -518,29 +616,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) @@ -568,9 +664,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: @@ -599,6 +692,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, @@ -608,8 +707,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] @@ -617,7 +722,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) @@ -630,18 +735,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) @@ -649,14 +751,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 @@ -692,7 +792,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 @@ -727,11 +827,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. """ @@ -767,10 +864,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]) @@ -781,106 +875,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: @@ -892,24 +888,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 @@ -917,19 +908,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): @@ -939,9 +923,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, \ @@ -954,9 +936,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) @@ -980,8 +960,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 @@ -992,6 +970,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 @@ -1006,8 +985,11 @@ def load_CNMF(filename:str, n_processes=1, dview=None): dview: multiprocessing or ipyparallel object used to set up parallelization, default None ''' + new_obj = CNMF(n_processes) - if os.path.splitext(filename)[1].lower() in ('.hdf5', '.h5'): + file_extension = os.path.splitext(filename)[1].lower() + + if file_extension in ('.hdf5', '.h5'): filename = caiman.paths.fn_relocated(filename) for key, val in load_dict_from_hdf5(filename).items(): if key == 'params': @@ -1034,9 +1016,8 @@ def load_CNMF(filename:str, n_processes=1, dview=None): setattr(new_obj, key, val) if new_obj.estimates.dims is None or new_obj.estimates.dims == b'NoneType': new_obj.estimates.dims = new_obj.dims - elif os.path.splitext(filename)[1].lower() == '.nwb': - from pynwb import NWBHDF5IO - with NWBHDF5IO(filename, 'r') as io: + elif file_extension == '.nwb': + with pynwb.NWBHDF5IO(filename, 'r') as io: nwb = io.read() ophys = nwb.processing['ophys'] rrs_group = ophys.data_interfaces['Fluorescence'].roi_response_series @@ -1058,6 +1039,7 @@ def load_CNMF(filename:str, n_processes=1, dview=None): else: b = None #np.zeros(mov.shape[1:]) f = None + estims = Estimates(A=A, b=b, C=C, f=f) estims.YrA = ophys.data_interfaces['residuals'].data[:].T @@ -1099,6 +1081,6 @@ def load_CNMF(filename:str, n_processes=1, dview=None): setattr(new_obj, 'estimates', estims) else: - raise NotImplementedError('unsupported file extension') + raise NotImplementedError(f'Unsupported file extension {file_extension}') return new_obj diff --git a/caiman/source_extraction/cnmf/deconvolution.py b/caiman/source_extraction/cnmf/deconvolution.py index 3322dd8d9..fd4d684b8 100644 --- a/caiman/source_extraction/cnmf/deconvolution.py +++ b/caiman/source_extraction/cnmf/deconvolution.py @@ -1001,8 +1001,8 @@ def estimate_time_constant(fluor, p=2, sn=None, lags=5, fudge_factor=1.): xc = axcov(fluor, lags) xc = xc[:, np.newaxis] - A = scipy.linalg.toeplitz(xc[lags + np.arange(lags)], - xc[lags + np.arange(p)]) - sn**2 * np.eye(lags, p) + A = scipy.linalg.toeplitz(c=np.ravel(xc[lags + np.arange(lags)]), + r=np.ravel(xc[lags + np.arange(p)])) - sn**2 * np.eye(lags, p) g = np.linalg.lstsq(A, xc[lags + 1:], rcond=None)[0] gr = np.roots(np.concatenate([np.array([1]), -g.flatten()])) gr = (gr + gr.conjugate()) / 2. @@ -1041,14 +1041,15 @@ def GetSn(fluor, range_ff=[0.25, 0.5], method='logmexp'): ind2 = ff < range_ff[1] ind = np.logical_and(ind1, ind2) Pxx_ind = Pxx[ind] - sn = { - 'mean': lambda Pxx_ind: np.sqrt(np.mean(Pxx_ind / 2)), - 'median': lambda Pxx_ind: np.sqrt(np.median(Pxx_ind / 2)), - 'logmexp': lambda Pxx_ind: np.sqrt(np.exp(np.mean(np.log(Pxx_ind / 2)))) - }[method](Pxx_ind) - - return sn + if method == 'mean': + return np.sqrt(np.mean(Pxx_ind / 2)) + elif method == 'median': + return np.sqrt(np.median(Pxx_ind / 2)) + elif method == 'logmexp': + return np.sqrt(np.exp(np.mean(np.log(Pxx_ind / 2)))) + else: + raise Exception('Invalid method requested for GetSn') def axcov(data, maxlag=5): """ 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..b63e3331e 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -14,8 +14,6 @@ from multiprocessing import current_process import numpy as np import scipy -import scipy.ndimage as nd -from scipy.ndimage import center_of_mass, correlate import scipy.sparse as spr from skimage.morphology import disk from sklearn.decomposition import NMF, FastICA @@ -148,7 +146,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 +167,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 +264,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. @@ -282,7 +293,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter if gSiz is None: gSiz = 2 * (np.asarray(gSig) + .5).astype(int) + 1 - d, T = np.shape(Y)[:-1], np.shape(Y)[-1] + d, T = Y.shape[:-1], Y.shape[-1] # rescale according to downsampling factor gSig = np.asarray(gSig, dtype=float) / ssub gSiz = np.round(np.asarray(gSiz) / ssub).astype(int) @@ -318,11 +329,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,16 +363,15 @@ 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] + K = Ain.shape[-1] if Ain.size > 0 and not center_psf and ssub != 1: @@ -397,7 +407,7 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter if Ain.size > 0: Cin = resize(Cin, [K, T]) center = np.asarray( - [center_of_mass(a.reshape(d, order='F')) for a in Ain.T]) + [scipy.ndimage.center_of_mass(a.reshape(d, order='F')) for a in Ain.T]) else: Cin = np.empty((K, T), dtype=np.float32) center = [] @@ -431,6 +441,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) @@ -442,7 +454,7 @@ def ICA_PCA(Y_ds, nr, sigma_smooth=(.5, .5, .5), truncate=2, fun='logcosh', m1 = m pca_comp = nr - T, d1, d2 = np.shape(m1) + T, d1, d2 = m1.shape d = d1 * d2 yr = np.reshape(m1, [T, d], order='F') @@ -583,30 +595,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 +689,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 +715,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 +727,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 +749,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,8 +763,8 @@ 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") - d = np.shape(Y) + logger.info(f"Greedy initialization of spatial and temporal components using spatial Gaussian filtering. Params={nmf_overrides}") + d = Y.shape Y[np.isnan(Y)] = 0 med = np.median(Y, axis=-1) Y = Y - med[..., np.newaxis] @@ -841,6 +855,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 +941,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) @@ -1035,12 +1061,10 @@ def imblur(Y, sig=5, siz=11, nDimBlur=None, kernel=None, opencv=True): h /= np.sqrt(h.dot(h)) shape = [1] * len(Y.shape) shape[i] = -1 - X = correlate(X, h.reshape(shape), mode='constant') + X = scipy.ndimage.correlate(X, h.reshape(shape), mode='constant') else: - X = correlate(Y, kernel[..., np.newaxis], mode='constant') - # for t in range(np.shape(Y)[-1]): - # X[:,:,t] = correlate(Y[:,:,t],kernel,mode='constant', cval=0.0) + X = scipy.ndimage.correlate(Y, kernel[..., np.newaxis], mode='constant') return X @@ -1081,13 +1105,13 @@ def hals(Y, A, C, b, f, bSiz=3, maxIter=5): """ # smooth the components - dims, T = np.shape(Y)[:-1], np.shape(Y)[-1] + dims, T = Y.shape[:-1], Y.shape[-1] K = A.shape[1] # number of neurons nb = b.shape[1] # number of background components if bSiz is not None: if isinstance(bSiz, (int, float)): bSiz = [bSiz] * len(dims) - ind_A = nd.filters.uniform_filter(np.reshape(A, + ind_A = scipy.ndimage.uniform_filter(np.reshape(A, dims + (K,), order='F'), size=bSiz + [0]) ind_A = np.reshape(ind_A > 1e-10, (np.prod(dims), K), order='F') else: @@ -1253,7 +1277,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 +1326,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 +1357,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..929b7e37b 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, @@ -252,13 +252,13 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, for jj, fff in enumerate(file_res): if fff is not None: idx_, shapes, A, b, C, f, S, bl, c1, neurons_sn, g, sn, _, YrA = fff - for _ in range(np.shape(b)[-1]): + for _ in range(b.shape[-1]): count_bgr += 1 A = A.tocsc() if del_duplicates: keep = [] - for ii in range(np.shape(A)[-1]): + for ii in range(A.shape[-1]): neuron_center = (np.array(scipy.ndimage.center_of_mass( A[:, ii].toarray().reshape(shapes, order='F'))) - np.array(shapes) / 2. + np.array(patch_centers[jj])) @@ -276,7 +276,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, file_res[jj][10] = g[keep] file_res[jj][-1] = YrA[keep] - # for ii in range(np.shape(A)[-1]): + # for ii in range(A.shape[-1]): # new_comp = A[:, ii] / np.sqrt(A[:, ii].power(2).sum()) # if new_comp.sum() > 0: # count += 1 @@ -332,7 +332,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, idx_ptr_B += list(b.indptr[1:] - b.indptr[:-1]) idx_tot_B.append(idx_[b.indices]) else: - for ii in range(np.shape(b)[-1]): + for ii in range(b.shape[-1]): b_tot.append(b[:, ii]) idx_tot_B.append(idx_) idx_ptr_B.append(len(idx_)) @@ -343,7 +343,7 @@ def run_CNMF_patches(file_name, shape, params, gnb=1, dview=None, else: # full background per patch F_tot = np.concatenate([F_tot, f]) - for ii in range(np.shape(A)[-1]): + for ii in range(A.shape[-1]): new_comp = A[:, ii] # / np.sqrt(A[:, ii].power(2).sum()) if new_comp.sum() > 0: a_tot.append(new_comp.toarray().flatten()) diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py index a01244472..f5c2281cf 100644 --- a/caiman/source_extraction/cnmf/merging.py +++ b/caiman/source_extraction/cnmf/merging.py @@ -139,7 +139,7 @@ def merge_components(Y, A, b, C, R, f, S, sn_pix, temporal_params, if R is None: R = np.zeros_like(C) - [d, t] = np.shape(Y) + d, t = Y.shape # find graph of overlapping spatial components A_corr = scipy.sparse.triu(A.T * A) @@ -169,7 +169,7 @@ def merge_components(Y, A, b, C, R, f, S, sn_pix, temporal_params, list_conxcomp = np.asarray(list_conxcomp_initial).T if list_conxcomp.ndim > 1: - cor = np.zeros((np.shape(list_conxcomp)[1], 1)) + cor = np.zeros((list_conxcomp.shape[1], 1)) for i in range(np.size(cor)): fm = np.where(list_conxcomp[:, i])[0] for j1 in range(np.size(fm)): diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py index 55b2828ba..14829d356 100644 --- a/caiman/source_extraction/cnmf/online_cnmf.py +++ b/caiman/source_extraction/cnmf/online_cnmf.py @@ -25,6 +25,7 @@ from scipy.sparse import coo_matrix, csc_matrix, spdiags, hstack from scipy.stats import norm from sklearn.decomposition import NMF +from skimage.morphology import disk from sklearn.preprocessing import normalize import tensorflow as tf from time import time @@ -47,9 +48,9 @@ from caiman.source_extraction.cnmf.utilities import (update_order, peak_local_max, decimation_matrix, gaussian_filter, uniform_filter) import caiman.summary_images -from caiman.utils.utils import save_dict_to_hdf5, load_dict_from_hdf5, parmap, load_graph -from caiman.utils.stats import pd_solve from caiman.utils.nn_models import (fit_NL_model, create_LN_model, quantile_loss, rate_scheduler) +from caiman.utils.stats import pd_solve +from caiman.utils.utils import save_dict_to_hdf5, load_dict_from_hdf5, parmap, load_graph try: cv2.setNumThreads(0) @@ -62,6 +63,10 @@ except: def profile(a): return a +#TODO If we ever get a chance, it would make sense to refactor CNMF and OnACID to have a +# parent class, as OnACID started as a copy of the CNMF codebase and they have similar +# APIs and intent, just a very different execution strategy. It would take some thought +# on how to do this without breaking compatibility, and on how big the parent class might be. class OnACID(object): """ Source extraction of streaming data using online matrix factorization. @@ -115,6 +120,37 @@ 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). + + def __str__(self): + ret = f"Caiman OnACID 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 OnACID 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 @profile def _prepare_object(self, Yr, T, new_dims=None, idx_components=None): @@ -149,7 +185,7 @@ def _prepare_object(self, Yr, T, new_dims=None, idx_components=None): if Yr.shape[-1] != self.params.get('online', 'init_batch'): raise Exception( - 'The movie size used for initialization does not match with the minibatch size') + 'The movie size used for initialization does not match the minibatch size') if new_dims is not None: @@ -303,6 +339,7 @@ def _prepare_object(self, Yr, T, new_dims=None, idx_components=None): self.estimates.rho_buf = RingBuffer(self.estimates.rho_buf, self.params.get('online', 'minibatch_shape')) self.estimates.sv = np.sum(self.estimates.rho_buf.get_last_frames( min(self.params.get('online', 'init_batch'), self.params.get('online', 'minibatch_shape')) - 1), 0) + self.estimates.AtA = (self.estimates.Ab.T.dot(self.estimates.Ab)).toarray() self.estimates.AtY_buf = self.estimates.Ab.T.dot(self.estimates.Yr_buf.T) self.estimates.groups = list(map(list, update_order(self.estimates.Ab)[0])) @@ -351,7 +388,6 @@ def _prepare_object(self, Yr, T, new_dims=None, idx_components=None): self.loaded_model = loaded_model if self.is1p: - from skimage.morphology import disk radius = int(round(self.params.get('init', 'ring_size_factor') * self.params.get('init', 'gSiz')[0] / float(ssub_B))) ring = disk(radius + 1) @@ -408,7 +444,7 @@ def get_indices_of_pixels_on_ring(self, pixel): @profile def fit_next(self, t, frame_in, num_iters_hals=3): """ - This method fits the next frame using the CaImAn online algorithm and + This method fits the next frame using the online algorithm and updates the object. Does NOT perform motion correction, see ``mc_next()`` Args @@ -421,6 +457,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 +527,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 +559,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, @@ -1041,6 +1078,7 @@ def mc_next(self, t: int, frame: np.ndarray) -> np.ndarray: templ += B else: templ = self.estimates.Ab.dot(self.estimates.C_on[:self.M, t-1]) + templ = templ.reshape(self.params.get('data', 'dims'), order='F') if self.params.get('online', 'normalize'): templ *= self.img_norm @@ -1084,8 +1122,8 @@ def mc_next(self, t: int, frame: np.ndarray) -> np.ndarray: def fit_online(self, **kwargs): """Implements the caiman online algorithm on the list of files fls. The - files are taken in alpha numerical order and are assumed to each have - the same number of frames (except the last one that can be shorter). + files are read in alphanumerical order and are assumed to each have + the same number of frames (except for the last, which can be shorter). Caiman online is initialized using the seeded or bare initialization methods. @@ -1106,9 +1144,8 @@ def fit_online(self, **kwargs): additional parameters used to modify self.params.online'] see options.['online'] for details - Returns: - self (results of caiman online) """ + logger = logging.getLogger("caiman") self.t_init = -time() fls = self.params.get('data', 'fnames') @@ -1148,6 +1185,7 @@ def fit_online(self, **kwargs): self.params.set('ring_CNN', {'path_to_model': path_to_model}) else: model_LN = None + epochs = self.params.get('online', 'epochs') self.initialize_online(model_LN=model_LN) self.t_init += time() @@ -1294,13 +1332,12 @@ def fit_online(self, **kwargs): self.estimates.C_on = self.estimates.C_on[:self.M] self.estimates.noisyC = self.estimates.noisyC[:self.M] - return self - def create_frame(self, frame_cor, show_residuals=True, resize_fact=3, transpose=True): if show_residuals: caption = 'Corr*PSNR buffer' if self.params.get('online', 'use_corr_img') else 'Mean Residual Buffer' else: caption = 'Identified Components' + captions = ['Raw Data', 'Inferred Activity', caption, 'Denoised Data'] self.dims = self.estimates.dims self.captions = captions @@ -1323,6 +1360,7 @@ def create_frame(self, frame_cor, show_residuals=True, resize_fact=3, transpose= self.estimates.W.dot(bc2))).reshape(self.dims, order='F') else: bgkrnd_frame = b.dot(f[:, self.t - 1]).reshape(self.dims, order='F') # denoised frame (components + background) + denoised_frame = comps_frame + bgkrnd_frame denoised_frame = (denoised_frame.copy() - self.bnd_Y[0])/np.diff(self.bnd_Y) comps_frame = (comps_frame.copy() - self.bnd_AC[0])/np.diff(self.bnd_AC) @@ -1340,7 +1378,7 @@ def create_frame(self, frame_cor, show_residuals=True, resize_fact=3, transpose= else: all_comps = np.array(A.sum(-1)).reshape(self.dims, order='F') fac = 2 - #all_comps = (all_comps.copy() - self.bnd_Y[0])/np.diff(self.bnd_Y) + all_comps = np.minimum(np.maximum(all_comps, 0)*fac, 1) # spatial shapes frame_comp_1 = cv2.resize(np.concatenate([frame_plot, all_comps * 1.], axis=-1), @@ -1774,7 +1812,7 @@ def refit(o, c): def init_shapes_and_sufficient_stats(Y, A, C, b, f, W=None, b0=None, ssub_B=1, bSiz=3, downscale_matrix=None, upscale_matrix=None): # smooth the components - dims, T = np.shape(Y)[:-1], np.shape(Y)[-1] + dims, T = Y.shape[:-1], Y.shape[-1] K = A.shape[1] # number of neurons if W is None: nb = b.shape[1] # number of background components @@ -2147,7 +2185,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, gHalf = np.array(gSiz) // 2 # number of total components (including background) - M = np.shape(Ab)[-1] + M = Ab.shape[-1] N = M - gnb # number of components (without background) if corr_img is None: @@ -2185,7 +2223,7 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, for ij in ijSig]), dims, order='F').ravel() cin_circ = cin.get_ordered() - useOASIS = False # whether to use faster OASIS for cell detection + useOASIS = False # whether to use faster OASIS for cell detection FIXME don't hardcode things internally like this accepted = True # flag indicating new component has not been rejected yet if Ab_dense is None: @@ -2241,18 +2279,32 @@ def update_num_components(t, sv, Ab, Cf, Yres_buf, Y_buf, rho_buf, ind_new.append(ijSig) if oases is not None: - if not useOASIS: + if not useOASIS: # FIXME bad variable name, also hardcoded? # lambda from Selesnick's 3*sigma*|K| rule # use noise estimate from init batch or use std_rr? # sn_ = sqrt((ain**2).dot(sn[indices]**2)) / sqrt(1 - g**2) + # The one-liner below was too hard to read -- breaking it apart for legibility + # (and also to make it easier to temporarily add assertions with np.isscalar and + # unpack size-1 arrays that are no longer ok in newer versions of numpy/scipy) sn_ = std_rr - oas = OASIS(np.ravel(g)[0], 3 * sn_ / - (sqrt(1 - g**2) if np.size(g) == 1 else - sqrt((1 + g[1]) * ((1 - g[1])**2 - g[0]**2) / (1 - g[1]))) - if s_min == 0 else 0, - s_min, num_empty_samples=t + - 1 - len(cin_res), - g2=0 if np.size(g) == 1 else g[1]) + if np.size(sn_) == 1: + sn_ = np.ravel(std_rr)[0] + else: + logger.warning("std_rr has more dimensionality than expected and this may lead to problems") + sn_ = std_rr + + oasis_g = np.ravel(g)[0] + if s_min != 0: + oasis_lambda = 0 + elif np.size(g) == 1: + oasis_lambda = 3 * sn_ / (sqrt(1 - np.ravel(g)[0]**2)) + else: + oasis_lambda = 3 * sn_ / sqrt((1 + np.ravel(g)[1]) * ((1 - np.ravel(g)[1])**2 - np.ravel(g)[0]**2) / (1 - np.ravel(g)[1])) + + oasis_ne = t + 1 - len(cin_res) + oasis_g2 = 0 if np.size(g) == 1 else np.ravel(g)[1] + + oas = OASIS(oasis_g, oasis_lambda, s_min, num_empty_samples=oasis_ne, g2=oasis_g2) for yt in cin_res: oas.fit_next(yt) @@ -2440,7 +2492,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 +2527,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..4c95e76ae 100644 --- a/caiman/source_extraction/cnmf/params.py +++ b/caiman/source_extraction/cnmf/params.py @@ -1,10 +1,10 @@ #!/usr/bin/env python +import importlib.metadata import json import logging import numpy as np import os -import pkg_resources from pprint import pformat import scipy from scipy.ndimage import generate_binary_structure, iterate_structure @@ -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 @@ -445,7 +453,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), If set to False, a list of submatrices is saved (typically faster). init_batch: int, default: 200, - length of mini batch used for initialization + length of mini batch used for initialization (must not exceed frame count on first file) init_method: 'bare'|'cnmf'|'seeded', default: 'bare', initialization method @@ -664,7 +672,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1), 'decay_time': decay_time, 'dxy': dxy, 'var_name_hdf5': var_name_hdf5, - 'caiman_version': pkg_resources.get_distribution('caiman').version, + 'caiman_version': importlib.metadata.version('caiman'), 'last_commit': None } @@ -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/pre_processing.py b/caiman/source_extraction/cnmf/pre_processing.py index d7eef9609..a0bc64b4d 100644 --- a/caiman/source_extraction/cnmf/pre_processing.py +++ b/caiman/source_extraction/cnmf/pre_processing.py @@ -114,7 +114,7 @@ def get_noise_welch(Y, noise_range=[0.25, 0.5], noise_method='logmexp', Y[..., int(T // 2 - max_num_samples_fft / 3 / 2): int(T // 2 + max_num_samples_fft / 3 / 2)], Y[..., -max_num_samples_fft // 3:]), axis=-1) - T = np.shape(Y)[-1] + T = Y.shape[-1] ff, Pxx = scipy.signal.welch(Y) Pxx = Pxx[..., (ff >= noise_range[0]) & (ff <= noise_range[1])] sn = { @@ -156,7 +156,7 @@ def get_noise_fft(Y, noise_range=[0.25, 0.5], noise_method='logmexp', max_num_sa Y[..., int(T // 2 - max_num_samples_fft / 3 / 2) :int(T // 2 + max_num_samples_fft / 3 / 2)], Y[..., -max_num_samples_fft // 3:]), axis=-1) - T = np.shape(Y)[-1] + T = Y.shape[-1] # we create a map of what is the noise on the FFT space ff = np.arange(0, 0.5 + 1. / T, 1. / T) @@ -396,7 +396,7 @@ def estimate_time_constant(Y, sn, p=None, lags=5, include_noise=False, pixels=No if p is None: raise Exception("You need to define p") if pixels is None: - pixels = np.arange(np.size(Y) // np.shape(Y)[-1]) + pixels = np.arange(np.size(Y) // Y.shape[-1]) npx = len(pixels) lags += p diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index 0ba5f2b4e..6e290b753 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -36,8 +36,8 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, min_size=3, max_size=8, dist=3, normalize_yyt_one=True, method_exp='dilate', expandCore=None, dview=None, n_pixels_per_process=128, - medw=(3, 3), thr_method='max', maxthr=0.1, - nrgthr=0.9999, extract_cc=True, b_in=None, + medw=(3, 3), thr_method='max', maxthr=0.1, nrgthr=0.9999, extract_cc=True, + b_in=None, se=np.ones((3, 3), dtype=int), ss=np.ones((3, 3), dtype=int), nb=1, method_ls='lasso_lars', update_background_components=True, @@ -66,41 +66,36 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, spatial profile of background activity. If A_in is boolean then it defines the spatial support of A. Otherwise it is used to determine it through determine_search_location - b_in: np.ndarray - you can pass background as input, especially in the case of one background per patch, since it will update using hals + sn: [optional] float + noise associated with each pixel if known dims: [optional] tuple x, y[, z] movie dimensions - min_size: [optional] int - - max_size: [optional] int - - dist: [optional] int - - sn: [optional] float - noise associated with each pixel if known + min_size, max_size, dist: [optional] int + Parameters for computing_indicator() - backend [optional] str - 'ipyparallel', 'single_thread' - single_thread:no parallelization. It can be used with small datasets. - ipyparallel: uses ipython clusters and then send jobs to each of them + normalize_yyt_one: bool + whether to normalize the C and A matrices so that diag(C*C.T) are ones - n_pixels_per_process: [optional] int - number of pixels to be processed by each thread - - method: [optional] string + method_exp: [optional] string method used to expand the search for pixels 'ellipse' or 'dilate' expandCore: [optional] scipy.ndimage.morphology if method is dilate this represents the kernel used for expansion dview: view on ipyparallel client - you need to create an ipyparallel client and pass a view on the processors (client = Client(), dview=client[:]) + you need to create an ipyparallel client and pass a view on the processors (client = Client(), dview=client[:]) + + n_pixels_per_process: [optional] int + number of pixels to be processed by each thread medw, thr_method, maxthr, nrgthr, extract_cc, se, ss: [optional] Parameters for components post-processing. Refer to spatial.threshold_components for more details + b_in: np.ndarray + you can pass background as input, especially in the case of one background per patch, since it will update using hals + nb: [optional] int Number of background components @@ -109,9 +104,6 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, 'nnls_L0'. Nonnegative least square with L0 penalty 'lasso_lars' lasso lars function from scikit learn - normalize_yyt_one: bool - whether to normalize the C and A matrices so that diag(C*C.T) are ones - update_background_components:bool whether to update the background components in the spatial phase @@ -119,6 +111,8 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, whether to update the using a low rank approximation. In the False case all the nonzero elements of the background components are updated using hals (to be used with one background per patch) + num_blocks_per_run_spat:int + (Undocumented) Returns: A: np.ndarray @@ -133,24 +127,6 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, f: np.ndarray same as f_in except if empty component deleted. - Raises: - Exception 'You need to define the input dimensions' - - Exception 'Dimension of Matrix Y must be pixels x time' - - Exception 'Dimension of Matrix C must be neurons x time' - - Exception 'Dimension of Matrix f must be background comps x time ' - - Exception 'Either A or C need to be determined' - - Exception 'Dimension of Matrix A must be pixels x neurons' - - Exception 'You need to provide estimate of C and f' - - Exception 'Not implemented consistently' - - Exception "Failed to delete: " + folder """ logger = logging.getLogger("caiman") # TODO fix documentation on backend @@ -180,13 +156,13 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, A_in = caiman.utils.stats.csc_column_remove(A_in, list(ff)) C = np.delete(C, list(ff), 0) # update indices - ind_list = list(range(nr-np.size(ff))) + ind_list = list(range(nr - np.size(ff))) for i in ff: ind_list.insert(i, 0) ind_list = np.array(ind_list, dtype=int) ind2_ = [ind_list[np.setdiff1d(a,ff)] if len(a) else a for a in ind2_] - nr = np.shape(C)[0] + nr = C.shape[0] if normalize_yyt_one and C is not None: C = np.array(C) d_ = scipy.sparse.lil_matrix((nr, nr)) @@ -246,8 +222,11 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if any(ff < nr): A_ = caiman.utils.stats.csc_column_remove(A_, list(ff[ff < nr])) C = np.delete(C, list(ff[ff < nr]), 0) - ff -= nr + + nr_ = nr nr = nr - len(ff[ff < nr]) + ff -= nr_ # previous nr seem to be updated incorrectly + else: ff -= nr if update_background_components: @@ -271,7 +250,6 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if b_in is None: # update baseline based on residual - #b = np.fmax(Y_resf.dot(np.linalg.inv(f.dot(f.T))), 0) b = np.fmax(np.linalg.solve(f.dot(f.T), Y_resf.T), 0).T else: ind_b = [np.where(_b)[0] for _b in b_in.T] @@ -281,11 +259,7 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None, if b_in is None: raise Exception( 'If you set the update_background_components to True you have to pass them as input to update_spatial') - # try: - # b = np.delete(b_in, background_ff, 0) - # except NameError: b = b_in - # print(("--- %s seconds ---" % (time.time() - start_time))) logger.info('Updating done in ' + '{0}s'.format(str(time.time() - start_time).split(".")[0])) try: # clean up @@ -371,7 +345,7 @@ def regression_ipyparallel(pars): else: C = C_name - _, T = np.shape(C) # initialize values + _, T = C.shape # initialize values As = [] for y, px, idx_px_from_0 in zip(Y, idxs_Y, range(len(idxs_C))): c = C[idxs_C[idx_px_from_0], :] @@ -491,7 +465,7 @@ def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgth if ss is None: ss = np.ones((3,) * len(dims), dtype='uint8') # dims and nm of neurones - d, nr = np.shape(A) + d, nr = A.shape # instantiation of A thresh. #Ath = np.zeros((d, nr)) pars = [] @@ -783,17 +757,15 @@ def test(Y, A_in, C, f, n_pixels_per_process, nb): if len(A_in.shape) == 1: A_in = np.atleast_2d(A_in).T if A_in.shape[0] == 1: - raise Exception( - 'Dimension of Matrix A must be pixels x neurons ') + raise Exception('Dimension of Matrix A must be pixels x neurons ') - [d, T] = np.shape(Y) + [d, T] = Y.shape if A_in is None: - A_in = np.ones((d, np.shape(C)[1]), dtype=bool) + A_in = np.ones((d, C.shape[1]), dtype=bool) if n_pixels_per_process > d: - print('The number of pixels per process (n_pixels_per_process)' - ' is larger than the total number of pixels!! Decreasing suitably.') + print(f'The number of pixels per process (n_pixels_per_process:{n_pixels_per_process}) is larger than the total number of pixels. Decreasing suitably.') n_pixels_per_process = d if f is not None: @@ -850,7 +822,7 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8, d1, d2 = dims elif len(dims) == 3: d1, d2, d3 = dims - d, nr = np.shape(A) + d, nr = A.shape A = csc_matrix(A) # dist_indicator = scipy.sparse.lil_matrix((d, nr),dtype= np.float32) # dist_indicator = scipy.sparse.csc_matrix((d, nr), dtype=np.float32) @@ -970,9 +942,9 @@ def construct_dilate_parallel(pars): return dist_indicator_i def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size, dist, expandCore, dview): - """compute the indices of the distance from the cm to search for the spatial component (calling determine_search_location) + """Compute the indices of the distance from the cm to search for the spatial component (calling determine_search_location) - does this by following an ellipse from the cm or doing a step by step dilatation around the cm + Uses strategy of following an ellipse from the cm or doing a step by step dilatation around the cm if it doesn't follow the rules it will throw an exception that is not supposed to be caught by spatial. @@ -980,57 +952,49 @@ def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size, Y: np.ndarray (2D or 3D) movie, raw data in 2D or 3D (pixels x time). - C: np.ndarray - calcium activity of each neuron. - - f: np.ndarray - temporal profile of background activity. - A_in: np.ndarray spatial profile of background activity. If A_in is boolean then it defines the spatial support of A. Otherwise it is used to determine it through determine_search_location - n_pixels_per_process: [optional] int - number of pixels to be processed by each thread + C: np.ndarray + calcium activity of each neuron. - min_size: [optional] int + f: np.ndarray + temporal profile of background activity. - max_size: [optional] int + nb: int + Number of background components - dist: [optional] int + method: string + method used to expand the search for pixels 'ellipse' or 'dilate' dims: [optional] tuple x, y[, z] movie dimensions - method: [optional] string - method used to expand the search for pixels 'ellipse' or 'dilate' - - expandCore: [optional] scipy.ndimage.morphology - if method is dilate this represents the kernel used for expansion + min_size: int + max_size: int - Returns: - same: - but reshaped and tested + dist: int - Raises: - Exception 'You need to define the input dimensions' + expandCore: scipy.ndimage.morphology (or None) + A numerical kernel used for expansion, only used if method is dilate - Exception 'Dimension of Matrix Y must be pixels x time' + dview (parallelism handler object) - Exception 'Dimension of Matrix C must be neurons x time' - - Exception 'Dimension of Matrix f must be background comps x time' - - Exception 'Either A or C need to be determined' - - Exception 'Dimension of Matrix A must be pixels x neurons' - - Exception 'You need to provide estimate of C and f' - - Exception 'Not implemented consistently' - - Exception 'Failed to delete: " + folder' + Returns: + _ind2: Indices (not documented) + + nr: int + Count of neurons + C + Reshaped/tested version of C (calcium activity of each neuron) + f + Reshaped/tested version of f (temporal profile of background) + b + Either a reshaped/tested version of b, a freshly computed one if C is None, or in some cases None + A_in + Reshaped/tested version of A_in (spatial profile of background activity) """ if A_in.dtype == bool: @@ -1043,7 +1007,7 @@ def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size, px = (np.sum(dist_indicator, axis=1) > 0) not_px = ~px - if nb>1: + if nb > 1: f = NMF(nb, init='nndsvda').fit(np.maximum(Y[not_px, :], 0)).components_ else: if Y.shape[-1] < 30000: @@ -1062,7 +1026,7 @@ def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size, C = np.maximum(csr_matrix(dist_indicator_av.T).dot( Y) - dist_indicator_av.T.dot(b).dot(f), 0) A_in = scipy.sparse.coo_matrix(A_in.astype(np.float32)) - nr, _ = np.shape(C) # number of neurons + nr, _ = C.shape # number of neurons ind2_ = [np.hstack((np.where(iid_)[0], nr + np.arange(f.shape[0]))) if np.size(np.where(iid_)[0]) > 0 else [] for iid_ in dist_indicator] @@ -1070,7 +1034,7 @@ def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size, if C is None: raise Exception('You need to provide estimate of C and f') - nr, _ = np.shape(C) # number of neurons + nr = C.shape[0] # number of neurons if b is None: dist_indicator = determine_search_location( diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 1ee1ee65a..195744188 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -34,9 +34,9 @@ def make_G_matrix(T, g): if isinstance(g, np.ndarray): if len(g) == 1 and g < 0: g = 0 - gs = np.matrix(np.hstack((1, -(g[:]).T))) - ones_ = np.matrix(np.ones((T, 1))) - G = spdiags((ones_ * gs).T, list(range(0, -len(g) - 1, -1)), T, T) + gs = np.array([np.hstack((1, -(g[:]).T))]) + ones_ = np.ones((T, 1)) + G = spdiags((ones_ @ gs).T, list(range(0, -len(g) - 1, -1)), T, T) return G else: @@ -49,7 +49,7 @@ def constrained_foopsi_parallel(arg_in): """ Ytemp, nT, jj_, bl, c1, g, sn, argss = arg_in - T = np.shape(Ytemp)[0] + T = Ytemp.shape[0] cc_, cb_, c1_, gn_, sn_, sp_, lam_ = caiman.source_extraction.cnmf.deconvolution.constrained_foopsi( Ytemp, bl=bl, c1=c1, g=g, sn=sn, **argss) gd_ = np.max(np.real(np.roots(np.hstack((1, -gn_.T))))) @@ -174,8 +174,8 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N raise Exception("You have to provide a value for p") # INITIALIZATION OF VARS - d, T = np.shape(Y) - nr = np.shape(A)[-1] + d, T = Y.shape + nr = A.shape[-1] if b is not None: # if b.shape[0] < b.shape[1]: # b = b.T @@ -193,7 +193,7 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N sn = np.repeat(None, nr) A = scipy.sparse.hstack((A, b)).tocsc() - S = np.zeros(np.shape(Cin)) + S = np.zeros(Cin.shape) Cin = np.vstack((Cin, fin)) C = Cin.copy() nA = np.ravel(A.power(2).sum(axis=0)) + np.finfo(np.float32).eps @@ -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/cnmf/utilities.py b/caiman/source_extraction/cnmf/utilities.py index 1a935ab67..8126438d6 100644 --- a/caiman/source_extraction/cnmf/utilities.py +++ b/caiman/source_extraction/cnmf/utilities.py @@ -572,7 +572,7 @@ def fast_prct_filt(input_data, level=8, frames_window=1000): """ data = np.atleast_2d(input_data).copy() - T = np.shape(data)[-1] + T = data.shape[-1] downsampfact = frames_window elm_missing = int(np.ceil(T * 1.0 / downsampfact) @@ -580,7 +580,7 @@ def fast_prct_filt(input_data, level=8, frames_window=1000): padbefore = int(np.floor(elm_missing / 2.)) padafter = int(np.ceil(elm_missing / 2.)) tr_tmp = np.pad(data.T, ((padbefore, padafter), (0, 0)), mode='reflect') - numFramesNew, num_traces = np.shape(tr_tmp) + numFramesNew, num_traces = tr_tmp.shape # compute baseline quickly tr_BL = np.reshape(tr_tmp, (downsampfact, int(numFramesNew / downsampfact), @@ -743,8 +743,8 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_ else: A = np.array(A) - d1, d2 = np.shape(Cn) - d, nr = np.shape(A) + d1, d2 = Cn.shape + d, nr = A.shape if max_number is None: max_number = nr @@ -760,9 +760,9 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_ cumEn /= cumEn[-1] Bvec = np.zeros(d) Bvec[indx] = cumEn - Bmat[i] = np.reshape(Bvec, np.shape(Cn), order='F') + Bmat[i] = np.reshape(Bvec, Cn.shape, order='F') - T = np.shape(Y)[-1] + T = Y.shape[-1] plt.close() fig = plt.figure() @@ -792,7 +792,7 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_ y3_tiny = Y[coords_y[0]:coords_y[-1] + 1, coords_x[0]:coords_x[-1] + 1, :] - dy_sz, dx_sz = np.shape(a3_tiny)[:-1] + dy_sz, dx_sz = a3_tiny.shape[:-1] y2_tiny = np.reshape(y3_tiny, (dx_sz * dy_sz, T), order='F') a2_tiny = np.reshape(a3_tiny, (dx_sz * dy_sz, nr), order='F') y2_res = y2_tiny - a2_tiny.dot(C) @@ -812,7 +812,7 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_ cumEn /= cumEn[-1] Bvec = np.zeros(d) Bvec[indx] = cumEn - bmat = np.reshape(Bvec, np.shape(Cn), order='F') + bmat = np.reshape(Bvec, Cn.shape, order='F') plt.contour(y, x, bmat, [thr]) plt.pause(.01) @@ -869,7 +869,7 @@ def update_order(A, new_a=None, prev_list=None, method='greedy'): Written by Eftychios A. Pnevmatikakis, Simons Foundation, 2015 ''' - K = np.shape(A)[-1] + K = A.shape[-1] if new_a is None and prev_list is None: if method == 'greedy': @@ -926,9 +926,9 @@ def order_components(A, C): A = np.array(A.todense()) nA2 = np.sqrt(np.sum(A**2, axis=0)) K = len(nA2) - A = np.array(np.matrix(A) * spdiags(1./nA2, 0, K, K)) + A = np.array(A @ spdiags(1./nA2, 0, K, K)) nA4 = np.sum(A**4, axis=0)**0.25 - C = np.array(spdiags(nA2, 0, K, K) * np.matrix(C)) + C = np.array(spdiags(nA2, 0, K, K) @ C) mC = np.ndarray.max(np.array(C), axis=1) srt = np.argsort(nA4 * mC)[::-1] A_or = A[:, srt] * spdiags(nA2[srt], 0, K, K) @@ -941,7 +941,7 @@ def update_order_random(A, flag_AA=True): randomized partitions of non-overlapping components """ - K = np.shape(A)[-1] + K = A.shape[-1] if flag_AA: AA = A.copy() else: @@ -992,7 +992,7 @@ def update_order_greedy(A, flag_AA=True): Author: Eftychios A. Pnevmatikakis, Simons Foundation, 2017 """ - K = np.shape(A)[-1] + K = A.shape[-1] parllcomp:list = [] for i in range(K): new_list = True diff --git a/caiman/source_extraction/volpy/mrcnn/model.py b/caiman/source_extraction/volpy/mrcnn/model.py index 2dee5f23f..ac8b8dd54 100644 --- a/caiman/source_extraction/volpy/mrcnn/model.py +++ b/caiman/source_extraction/volpy/mrcnn/model.py @@ -26,8 +26,8 @@ from ..mrcnn import utils # Requires TensorFlow 2.0+ -from distutils.version import LooseVersion -assert LooseVersion(tf.__version__) >= LooseVersion("2.0") +from packaging.version import Version +assert Version(tf.__version__) >= Version("2.0") tf.compat.v1.disable_eager_execution() diff --git a/caiman/source_extraction/volpy/mrcnn/utils.py b/caiman/source_extraction/volpy/mrcnn/utils.py index 716ef3478..ac4e87746 100644 --- a/caiman/source_extraction/volpy/mrcnn/utils.py +++ b/caiman/source_extraction/volpy/mrcnn/utils.py @@ -7,21 +7,21 @@ Written by Waleed Abdulla """ -import sys -import os import logging import math -import random import numpy as np -import tensorflow as tf +import os +from packaging.version import Version +import random import scipy +import shutil import skimage.color import skimage.io import skimage.transform +import sys +import tensorflow as tf import urllib.request -import shutil import warnings -from distutils.version import LooseVersion # URL from which to download the latest COCO trained weights COCO_MODEL_URL = "https://github.com/matterport/Mask_RCNN/releases/download/v2.0/mask_rcnn_coco.h5" @@ -893,7 +893,7 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True, of skimage. This solves the problem by using different parameters per version. And it provides a central place to control resizing defaults. """ - if LooseVersion(skimage.__version__) >= LooseVersion("0.14"): + if Version(skimage.__version__) >= Version("0.14"): # New in 0.14: anti_aliasing. Default it to False for backward # compatibility with skimage 0.13. return skimage.transform.resize( 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/summary_images.py b/caiman/summary_images.py index de4e40792..dc00ee521 100644 --- a/caiman/summary_images.py +++ b/caiman/summary_images.py @@ -206,7 +206,7 @@ def local_correlations(Y, eight_neighbours: bool = True, swap_dim: bool = True, if swap_dim: Y = np.transpose(Y, tuple(np.hstack((Y.ndim - 1, list(range(Y.ndim))[:-1])))) - rho = np.zeros(np.shape(Y)[1:]) + rho = np.zeros(Y.shape[1:]) w_mov = (Y - np.mean(Y, axis=0)) / np.std(Y, axis=0) rho_h = np.mean(np.multiply(w_mov[:, :-1, :], w_mov[:, 1:, :]), axis=0) @@ -214,7 +214,7 @@ def local_correlations(Y, eight_neighbours: bool = True, swap_dim: bool = True, # yapf: disable if order_mean == 0: - rho = np.ones(np.shape(Y)[1:]) + rho = np.ones(Y.shape[1:]) rho_h = rho_h rho_w = rho_w rho[:-1, :] = rho[:-1, :] * rho_h @@ -232,7 +232,7 @@ def local_correlations(Y, eight_neighbours: bool = True, swap_dim: bool = True, rho[:, :, :-1] = rho[:, :, :-1] + rho_d rho[:, :, 1:] = rho[:, :, 1:] + rho_d - neighbors = 6 * np.ones(np.shape(Y)[1:]) + neighbors = 6 * np.ones(Y.shape[1:]) neighbors[0] = neighbors[0] - 1 neighbors[-1] = neighbors[-1] - 1 neighbors[:, 0] = neighbors[:, 0] - 1 @@ -258,7 +258,7 @@ def local_correlations(Y, eight_neighbours: bool = True, swap_dim: bool = True, rho[1:, :-1] = rho[1:, :-1] + rho_d1**(order_mean) rho[:-1, 1:] = rho[:-1, 1:] + rho_d2**(order_mean) - neighbors = 8 * np.ones(np.shape(Y)[1:3]) + neighbors = 8 * np.ones(Y.shape[1:3]) neighbors[0, :] = neighbors[0, :] - 3 neighbors[-1, :] = neighbors[-1, :] - 3 neighbors[:, 0] = neighbors[:, 0] - 3 @@ -268,7 +268,7 @@ def local_correlations(Y, eight_neighbours: bool = True, swap_dim: bool = True, neighbors[-1, 0] = neighbors[-1, 0] + 1 neighbors[0, -1] = neighbors[0, -1] + 1 else: - neighbors = 4 * np.ones(np.shape(Y)[1:3]) + neighbors = 4 * np.ones(Y.shape[1:3]) neighbors[0, :] = neighbors[0, :] - 1 neighbors[-1, :] = neighbors[-1, :] - 1 neighbors[:, 0] = neighbors[:, 0] - 1 diff --git a/caiman/tests/comparison/create_gt.py b/caiman/tests/comparison/create_gt.py index 1954f7c4f..a3ec6d8e2 100644 --- a/caiman/tests/comparison/create_gt.py +++ b/caiman/tests/comparison/create_gt.py @@ -147,7 +147,7 @@ def create(): m_orig = cm.load(fname) min_mov = m_orig[:400].min() comp = comparison.Comparison() - comp.dims = np.shape(m_orig)[1:] + comp.dims = m_orig.shape[1:] ################ RIG CORRECTION ################# t1 = time.time() @@ -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..04de626ee 100644 --- a/caiman/tests/comparison_general.py +++ b/caiman/tests/comparison_general.py @@ -150,7 +150,7 @@ def test_general(): m_orig = cm.load(fname) min_mov = m_orig[:400].min() comp = comparison.Comparison() - comp.dims = np.shape(m_orig)[1:] + comp.dims = m_orig.shape[1:] ################ RIG CORRECTION ################# t1 = time.time() @@ -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_import.py b/caiman/tests/test_import.py index e87ce15f4..213e6003b 100644 --- a/caiman/tests/test_import.py +++ b/caiman/tests/test_import.py @@ -2,4 +2,4 @@ def test_base(): - return True + pass 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_temporal.py b/caiman/tests/test_temporal.py index 5049bf6a4..c5b7a2b41 100644 --- a/caiman/tests/test_temporal.py +++ b/caiman/tests/test_temporal.py @@ -11,7 +11,7 @@ def test_make_G_matrix(): G = cnmf.temporal.make_G_matrix(T, g) G = G.todense() # yapf: disable - true_G = np.matrix( + true_G = np.array( [[1., 0., 0., 0., 0., 0.], [-1., 1., 0., 0., 0., 0.], [-2., -1., 1., 0., 0., 0.], diff --git a/caiman/tests/test_toydata.py b/caiman/tests/test_toydata.py index c796eaddd..3c63334da 100644 --- a/caiman/tests/test_toydata.py +++ b/caiman/tests/test_toydata.py @@ -3,13 +3,13 @@ import numpy.testing as npt import numpy as np import os -from scipy.ndimage.filters import gaussian_filter +import scipy.ndimage +from caiman import save_memmap, load_memmap +from caiman.paths import fn_relocated, generate_fname_tot import caiman.source_extraction.cnmf.params from caiman.source_extraction import cnmf as cnmf from caiman.utils.visualization import get_contours -from caiman.paths import fn_relocated, generate_fname_tot -from caiman import save_memmap, load_memmap TOYDATA_DIMS = { 2: (20, 30), @@ -34,8 +34,8 @@ def gen_data(D=3, noise=.5, T=300, framerate=30, firerate=2.): trueA[tuple(centers[i]) + (i,)] = 1. tmp = np.zeros(dims) tmp[tuple(d // 2 for d in dims)] = 1. - z = np.linalg.norm(gaussian_filter(tmp, sig).ravel()) - trueA = 10 * gaussian_filter(trueA, sig + (0,)) / z + z = np.linalg.norm(scipy.ndimage.gaussian_filter(tmp, sig).ravel()) + trueA = 10 * scipy.ndimage.gaussian_filter(trueA, sig + (0,)) / z Yr = bkgrd + noise * np.random.randn(*(np.prod(dims), T)) + \ trueA.reshape((-1, 4), order='F').dot(trueC) return Yr, trueC, trueS, trueA, centers, dims @@ -95,7 +95,6 @@ def get_params_dicts(D: int): } } - def pipeline(D, params_dict, name): #%% GENERATE GROUND TRUTH DATA Yr, trueC, trueS, trueA, centers, dims = gen_data(D) @@ -118,7 +117,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..ddb74eb38 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: @@ -35,7 +35,7 @@ def _check_keys(checkdict: dict) -> None: """ for key in checkdict: - if isinstance(checkdict[key], scipy.io.matlab.mio5_params.mat_struct): + if isinstance(checkdict[key], scipy.io.matlab.mat_struct): checkdict[key] = _todict(checkdict[key]) @@ -47,7 +47,7 @@ def _todict(matobj) -> dict: ret = {} for strg in matobj._fieldnames: elem = matobj.__dict__[strg] - if isinstance(elem, scipy.io.matlab.mio5_params.mat_struct): + if isinstance(elem, scipy.io.matlab.mat_struct): ret[strg] = _todict(elem) else: ret[strg] = elem @@ -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/caiman/utils/visualization.py b/caiman/utils/visualization.py index 9f9a6f2e1..20367cbc2 100644 --- a/caiman/utils/visualization.py +++ b/caiman/utils/visualization.py @@ -81,8 +81,7 @@ def view_patches(Yr, A, C, b, f, d1, d2, YrA=None, secs=1): A2.data **= 2 nA2 = np.sqrt(np.array(A2.sum(axis=0))).squeeze() if YrA is None: - Y_r = np.array(A.T * np.matrix(Yr) - (A.T * np.matrix(b[:, np.newaxis])) * np.matrix( - f[np.newaxis]) - (A.T.dot(A)) * np.matrix(C) + C) + Y_r = np.array(A.T @ Yr - (A.T @ b[:, np.newaxis]) @ f[np.newaxis] - (A.T.dot(A)) @ C + C) else: Y_r = YrA + C @@ -162,9 +161,8 @@ def nb_view_patches(Yr, A, C, b, f, d1, d2, YrA=None, image_neurons=None, thr=0. f = np.squeeze(f) if YrA is None: Y_r = np.array(spdiags(1 / nA2, 0, nr, nr) * - (A.T * np.matrix(Yr) - - (A.T * np.matrix(b[:, np.newaxis])) * np.matrix(f[np.newaxis]) - - A.T.dot(A) * np.matrix(C)) + C) + (A.T @ Yr - + (A.T @ b[:, np.newaxis]) @ f[np.newaxis] - A.T.dot(A) @ C) + C) else: Y_r = C + YrA @@ -321,9 +319,7 @@ def hv_view_patches(Yr, A, C, b, f, d1, d2, YrA=None, image_neurons=None, denois if YrA is None: Y_r = np.array( spdiags(1 / nA2, 0, nr, nr) * - (A.T * np.matrix(Yr) - - (A.T * np.matrix(b[:, np.newaxis])) * np.matrix(f[np.newaxis]) - - A.T.dot(A) * np.matrix(C)) + C) + (A.T @ Yr - (A.T @ b[:, np.newaxis]) @ f[np.newaxis] - A.T.dot(A) @ C) + C) else: Y_r = C + YrA if image_neurons is None: @@ -400,7 +396,7 @@ def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False, slice_dim: if 'csc_matrix' not in str(type(A)): A = csc_matrix(A) - d, nr = np.shape(A) + d, nr = A.shape coordinates = [] @@ -574,8 +570,8 @@ def nb_view_patches3d(Y_r, A, C, dims, image_type='mean', Yr=None, plt.close() K = np.max([[len(cor['coordinates']) for cor in cc] for cc in coors]) - cc1 = np.nan * np.zeros(np.shape(coors) + (K,)) - cc2 = np.nan * np.zeros(np.shape(coors) + (K,)) + cc1 = np.nan * np.zeros(coors.shape + (K,)) + cc2 = np.nan * np.zeros(coors.shape + (K,)) for i, cor in enumerate(coors[0]): cc1[0, i, :len(cor['coordinates']) ] = cor['coordinates'][:, 0] + offset1 @@ -852,7 +848,7 @@ def nb_plot_contour(image, A, d1, d2, thr=None, thr_method='max', maxthr=0.2, nr p.circle(center[:, 1], center[:, 0], size=10, color="black", fill_color=None, line_width=2, alpha=1) if coordinates is None: - coors = get_contours(A, np.shape(image), thr, thr_method) + coors = get_contours(A, image.shape, thr, thr_method) else: coors = coordinates cc1 = [np.clip(cor['coordinates'][1:-1, 0], 0, d2) for cor in coors] @@ -1103,7 +1099,7 @@ def plot_contours(A, Cn, thr=None, thr_method='max', maxthr=0.2, nrgthr=0.9, dis plt.imshow(Cn, interpolation=None, cmap=cmap, vmin=vmin, vmax=vmax) if coordinates is None: - coordinates = get_contours(A, np.shape(Cn), thr, thr_method, swap_dim) + coordinates = get_contours(A, Cn.shape, thr, thr_method, swap_dim) for c in coordinates: v = c['coordinates'] c['bbox'] = [np.floor(np.nanmin(v[:, 1])), np.ceil(np.nanmax(v[:, 1])), diff --git a/demos/general/demo_behavior.py b/demos/general/demo_behavior.py index 6d441bb06..169aaad45 100755 --- a/demos/general/demo_behavior.py +++ b/demos/general/demo_behavior.py @@ -90,7 +90,7 @@ def main(): spfl = spatial_filter spfl = cm.movie(spfl[None, :, :]).resize( 1 / resize_fact, 1 / resize_fact, 1).squeeze() - max_x, max_y = np.add((min_x, min_y), np.shape(spfl)) + max_x, max_y = np.add((min_x, min_y), spfl.shape) mask[min_x:max_x, min_y:max_y] = spfl mask[mask < np.nanpercentile(spfl, 70)] = np.nan 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..a85d3c751 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. @@ -121,8 +121,7 @@ def main(): moviehandle.play(fr=40, q_max=99.5, magnification=4) # press q to exit # MEMORY MAPPING - do_memory_mapping = True - if do_memory_mapping: + if not cfg.no_memmap: border_to_0 = 0 if mc.border_nan == 'copy' else mc.border_to_0 # you can include the boundaries of the FOV if you used the 'copy' option # during motion correction, although be careful about the components near @@ -167,7 +166,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 +220,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) @@ -268,6 +267,7 @@ def main(): def handle_args(): parser = argparse.ArgumentParser(description="Demonstrate voltage imaging pipeline") + parser.add_argument("--no_memmap", action="store_true", help="Avoid memory mapping") parser.add_argument("--keep_logs", action="store_true", help="Keep temporary logfiles") parser.add_argument("--no_play", action="store_true", help="Do not display results") parser.add_argument("--cluster_backend", default="multiprocessing", help="Specify multiprocessing, ipyparallel, or single to pick an engine") 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..c28e89556 100644 --- a/environment-minimal.yml +++ b/environment-minimal.yml @@ -11,7 +11,6 @@ dependencies: - ipywidgets - matplotlib - moviepy -- pytest - numpy <2.0.0,>=1.26 - numpydoc - opencv @@ -19,6 +18,8 @@ dependencies: - pims - psutil - pynwb +- pyside6 +- pytest - scikit-image >=0.19.0 - scikit-learn >=1.2 - scipy >= 1.10.1 diff --git a/environment.yml b/environment.yml index 014c3532c..a4dd90e86 100644 --- a/environment.yml +++ b/environment.yml @@ -11,13 +11,12 @@ dependencies: - ipykernel - ipython - ipyparallel +- ipywidgets - jupyter - jupyter_bokeh - matplotlib - moviepy - mypy -- pytest -- pytest-cov - numpy <2.0.0,>=1.26 - numpydoc - opencv @@ -27,6 +26,9 @@ dependencies: - psutil - pynwb - pyqtgraph +- pyside6 +- pytest +- pytest-cov - 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..8b1eb271d 100755 --- a/setup.py +++ b/setup.py @@ -1,8 +1,8 @@ #!/usr/bin/env python -from setuptools import setup, find_packages import numpy as np import os +from setuptools import setup, find_packages import sys from Cython.Build import cythonize from setuptools.extension import Extension @@ -29,12 +29,12 @@ # # We can access these by using sys.prefix as the base of the directory and constructing from there. # Note that if python's packaging standards ever change the install base of data_files to be under the -# package that made them, we can switch to using the pkg_resources API. +# package that made them, we can switch to using a different API. 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