Merge pull request flatironinstitute#618 from PBLab/dev

Updates to fit_file

.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

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

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))

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

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)