diff --git a/.gitignore b/.gitignore index 73ea7a621..db8d64aa3 100644 --- a/.gitignore +++ b/.gitignore @@ -17,7 +17,7 @@ Icon ._* *.log *.err -*.*.out +*.*.out # Files that might appear in the root of a volume .DocumentRevisions-V100 @@ -75,6 +75,7 @@ wheels/ *.egg-info/ .installed.cfg *.egg +.vscode/ # PyInstaller # Usually these files are written by a python script from a template diff --git a/caiman/mmapping.py b/caiman/mmapping.py index 5dff1db12..c052f2778 100644 --- a/caiman/mmapping.py +++ b/caiman/mmapping.py @@ -21,6 +21,7 @@ import sys import tifffile from typing import Any, Dict, List, Optional, Tuple, Union +import pathlib import caiman as cm from caiman.paths import memmap_frames_filename @@ -34,6 +35,7 @@ def prepare_shape(mytuple: Tuple) -> Tuple: return tuple(map(lambda x: np.uint64(x), mytuple)) + #%% def load_memmap(filename: str, mode: str = 'r') -> Tuple[Any, Tuple, int]: """ Load a memory mapped file created by the function save_memmap @@ -56,24 +58,23 @@ def load_memmap(filename: str, mode: str = 'r') -> Tuple[Any, Tuple, int]: Raises: - Exception "Unknown file extension" + ValueError "Unknown file extension" """ - if ('.mmap' in filename): - # Strip path components and use CAIMAN_DATA/example_movies - # TODO: Eventually get the code to save these in a different dir - file_to_load = filename - filename = os.path.split(filename)[-1] - fpart = filename.split('_')[1:-1] # The filename encodes the structure of the map - d1, d2, d3, T, order = int(fpart[-9]), int(fpart[-7]), int(fpart[-5]), int(fpart[-1]), fpart[-3] - Yr = np.memmap(file_to_load, mode=mode, shape=prepare_shape((d1 * d2 * d3, T)), dtype=np.float32, order=order) - if d3 == 1: - return (Yr, (d1, d2), T) - else: - return (Yr, (d1, d2, d3), T) - else: + if pathlib.Path(filename).suffix != '.mmap': logging.error("Unknown extension for file " + str(filename)) - raise Exception('Unknown file extension (should be .mmap)') + raise ValueError('Unknown file extension (should be .mmap)') + # Strip path components and use CAIMAN_DATA/example_movies + # TODO: Eventually get the code to save these in a different dir + file_to_load = filename + filename = os.path.split(filename)[-1] + fpart = filename.split('_')[1:-1] # The filename encodes the structure of the map + d1, d2, d3, T, order = int(fpart[-9]), int(fpart[-7]), int(fpart[-5]), int(fpart[-1]), fpart[-3] + Yr = np.memmap(file_to_load, mode=mode, shape=prepare_shape((d1 * d2 * d3, T)), dtype=np.float32, order=order) + if d3 == 1: + return (Yr, (d1, d2), T) + else: + return (Yr, (d1, d2, d3), T) #%% @@ -386,9 +387,9 @@ def save_memmap(filenames: List[str], slices: slice object or list of slice objects slice can be used to select portion of the movies in time and x,y - directions. For instance - slices = [slice(0,200),slice(0,100),slice(0,100)] will take - the first 200 frames and the 100 pixels along x and y dimensions. + directions. For instance + slices = [slice(0,200),slice(0,100),slice(0,100)] will take + the first 200 frames and the 100 pixels along x and y dimensions. Returns: fname_new: the name of the mapped file, the format is such that the name will contain the frame dimensions and the number of frames diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py index 84eb845fd..7ed8dd51c 100644 --- a/caiman/source_extraction/cnmf/cnmf.py +++ b/caiman/source_extraction/cnmf/cnmf.py @@ -32,6 +32,8 @@ import psutil import scipy import sys +import glob +import pathlib from .estimates import Estimates from .initialization import initialize_components, compute_W @@ -46,6 +48,8 @@ from ...components_evaluation import estimate_components_quality from ...motion_correction import MotionCorrect from ...utils.utils import save_dict_to_hdf5, load_dict_from_hdf5 +from caiman import summary_images +from caiman import cluster try: cv2.setNumThreads(0) @@ -302,24 +306,28 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p self.estimates = Estimates(A=Ain, C=Cin, b=b_in, f=f_in, dims=self.params.data['dims']) - def fit_file(self, motion_correct=False, indices=[slice(None)]*2): + def fit_file(self, motion_correct=False, indices=None, include_eval=False): """ This method packages the analysis pipeline (motion correction, memory - mapping, patch based CNMF processing) in a single method that can be - called on a specific (sequence of) file(s). It is assumed that the CNMF - object already contains a params object where the location of the files - and all the relevant parameters have been specified. The method does - not perform the quality evaluation step. Consult demo_pipeline for an - example. + 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 + where the location of the files and all the relevant parameters have + been specified. The method will perform the last step, i.e. component + evaluation, if the flag "include_eval" is set to `True`. Args: motion_correct (bool) flag for performing motion correction indices (list of slice objects) perform analysis only on a part of the FOV + include_eval (bool) + flag for performing component evaluation Returns: cnmf object with the current estimates """ + if indices is None: + indices = (slice(None), slice(None)) fnames = self.params.get('data', 'fnames') if os.path.exists(fnames[0]): _, extension = os.path.splitext(fnames[0])[:2] @@ -328,6 +336,7 @@ def fit_file(self, motion_correct=False, indices=[slice(None)]*2): logging.warning("Error: File not found, with file list:\n" + fnames[0]) raise Exception('File not found!') + base_name = pathlib.Path(fnames[0]).stem + "_memmap_" if extension == '.mmap': fname_new = fnames[0] Yr, dims, T = mmapping.load_memmap(fnames[0]) @@ -345,17 +354,45 @@ def fit_file(self, motion_correct=False, indices=[slice(None)]*2): else: b0 = np.ceil(np.max(np.abs(mc.shifts_rig))).astype(np.int) self.estimates.shifts = mc.shifts_rig - b0 = 0 if self.params.get('motion', 'border_nan') is 'copy' else 0 - fname_new = mmapping.save_memmap(fname_mc, base_name='memmap_', order='C', + # TODO - b0 is currently direction inspecific, which can cause + # sub-optimal behavior. See + # https://github.com/flatironinstitute/CaImAn/pull/618#discussion_r313960370 + # for further details. + # b0 = 0 if self.params.get('motion', 'border_nan') is 'copy' else 0 + b0 = 0 + fname_new = mmapping.save_memmap(fname_mc, base_name=base_name, order='C', border_to_0=b0) else: - fname_new = mmapping.save_memmap(fnames, base_name='memmap_', order='C') - # now load the file + fname_new = mmapping.save_memmap(fnames, base_name=base_name, order='C') Yr, dims, T = mmapping.load_memmap(fname_new) images = np.reshape(Yr.T, [T] + list(dims), order='F') self.mmap_file = fname_new - return self.fit(images, indices=indices) + if not include_eval: + return self.fit(images, indices=indices) + + fit_cnm = self.fit(images, indices=indices) + Cn = summary_images.local_correlations(images, swap_dim=False) + Cn[np.isnan(Cn)] = 0 + fit_cnm.save(fname_new[:-5]+'_init.hdf5') + fit_cnm.params.change_params({'p': self.params.get('preprocess', 'p')}) + # %% RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution + cnm2 = fit_cnm.refit(images, dview=self.dview) + cnm2.estimates.evaluate_components(images, cnm2.params, dview=self.dview) + #%% update object with selected components + cnm2.estimates.select_components(use_object=True) + #%% Extract DF/F values + cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250) + cnm2.estimates.Cn = Cn + cnm2.save(cnm2.mmap_file[:-4] + 'hdf5') + + cluster.stop_server(dview=self.dview) + log_files = glob.glob('*_LOG_*') + for log_file in log_files: + os.remove(log_file) + + return cnm2 + def refit(self, images, dview=None): """ @@ -378,7 +415,7 @@ def refit(self, images, dview=None): cnm.mmap_file = self.mmap_file return cnm.fit(images) - def fit(self, images, indices=[slice(None), slice(None)]): + def fit(self, images, indices=(slice(None), slice(None))): """ This method uses the cnmf algorithm to find sources in data. it is calling every function from the cnmf folder @@ -404,6 +441,8 @@ 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)) diff --git a/caiman/source_extraction/cnmf/params.py b/caiman/source_extraction/cnmf/params.py index 55a3323a9..1f1e9488a 100644 --- a/caiman/source_extraction/cnmf/params.py +++ b/caiman/source_extraction/cnmf/params.py @@ -572,7 +572,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': '1.5.2', + 'caiman_version': '1.5.3', 'last_commit': None, 'mmap_F': None, 'mmap_C': None diff --git a/caiman/tests/test_memmapping.py b/caiman/tests/test_memmapping.py new file mode 100644 index 000000000..6a3e6df9f --- /dev/null +++ b/caiman/tests/test_memmapping.py @@ -0,0 +1,75 @@ +import pathlib + +import numpy as np +import nose + +from caiman import mmapping +from caiman.paths import caiman_datadir + + +TWO_D_FNAME = ( + pathlib.Path(caiman_datadir()) + / "testdata/memmap__d1_10_d2_11_d3_1_order_F_frames_12_.mmap" +) +THREE_D_FNAME = ( + pathlib.Path(caiman_datadir()) + / "testdata/memmap__d1_10_d2_11_d3_13_order_F_frames_12_.mmap" +) + + +def test_load_raises_wrong_ext(): + fname = "a.mmapp" + try: + mmapping.load_memmap(fname) + except ValueError: + assert True + else: + assert False + + +def test_load_raises_multiple_ext(): + fname = "a.mmap.mma" + try: + mmapping.load_memmap(fname) + except ValueError: + assert True + else: + assert False + + +def setup_2d_mmap(): + np.memmap( + TWO_D_FNAME, mode="w+", dtype=np.float32, shape=(12, 10, 11, 13), order="F" + ) + + +def teardown_2d_mmap(): + TWO_D_FNAME.unlink() + + +def setup_3d_mmap(): + np.memmap( + THREE_D_FNAME, mode="w+", dtype=np.float32, shape=(12, 10, 11, 13), order="F" + ) + + +def teardown_3d_mmap(): + THREE_D_FNAME.unlink() + + +@nose.with_setup(setup_2d_mmap, teardown_2d_mmap) +def test_load_successful_2d(): + fname = pathlib.Path(caiman_datadir()) / "testdata" / TWO_D_FNAME + Yr, (d1, d2), T = mmapping.load_memmap(str(fname)) + assert (d1, d2) == (10, 11) + assert T == 12 + assert isinstance(Yr, np.memmap) + + +@nose.with_setup(setup_3d_mmap, teardown_3d_mmap) +def test_load_successful_3d(): + fname = pathlib.Path(caiman_datadir()) / "testdata" / THREE_D_FNAME + Yr, (d1, d2, d3), T = mmapping.load_memmap(str(fname)) + assert (d1, d2, d3) == (10, 11, 13) + assert T == 12 + assert isinstance(Yr, np.memmap)