+
+CaImAn comes with an experimental visual interface. To see an example on how use it, first load and run either
+* demo_OnACID_mesoscope.py
+* demo_caiman_basic.py
+* demo_pipeline.py
+
+Each of these demos will save a results file. You can then start the visual interface by running the following
+command from the base caiman folder (make sure you are within your caiman environment):
+```
+ipython caiman/gui/gui_pyqtgraph_layout.py
+```
+You will then be asked to load the file that is generated at the end of each file (ending in .hdf5)
+The associated .mmap files also need to be present in the same directory as the .hdf5
+
+A visual interface will appear, where you will be able to:
+* regulate gain and contrast for the background image (correlation image)
+* regulate the threshold over the spatial masks to visualize components contours
+* click on neurons and see the corresponding trace and mask
+* select subset of neurons based on different quality metrics
+* save the resulting selection to a file in hdf5 format
+
+More features will be added in the future.
diff --git a/README.md b/README.md
index 1f06abde3..afcee1d2f 100755
--- a/README.md
+++ b/README.md
@@ -34,6 +34,15 @@ A paper explaining most of the implementation details and benchmarking can be fo
All the results and figures of the paper can be regenerated using this package. For more information visit this [page](https://github.com/flatironinstitute/CaImAn/tree/master/use_cases/eLife_scripts).
+## New: Exporting results, GUI and NWB support (July 2019)
+
+You can now use the `save` method included in both the `CNMF` and `OnACID` objects to export the results (and parameters used) of your analysis. The results are saved in an HDF5 file that you can then load in a graphical user interface for more inspection. The GUI will allow you to inspect the results and modify the selected components based on the various quality metrics. For more information click [here](GUI.md)
+
+The [Neurodata Without Borders (NWB)](https://www.nwb.org/) file format is now supported by CaImAn. You read and analyze NWB files and can save the results of the analysis (`Estimates` object) back to the original NWB file. Consult this [demo](use_cases/NWB/demo_pipeline_NWB.py) for an example on how to use this feature.
+
+**To use CaImAn with these additional features you'll need to create a new environment following the usual instructions.**
+
+
## New: Removing Keras dependency (June 2019)
To circumvent a problem arising in Windows installation, we recently removed keras from the list of dependencies. Keras was being used to deploy the pretrained neural network models for component screening. The neural network models are being deployed through tensorflow and for that we included tensorflow compatible versions of the models inside the folder `model`. Existing users who already have keras in their environment will continue using keras as it is slightly faster. However if you create an environment without it you may want to either reinstall caimanmanager or simply copy the files `model/*.pb` into the folder `caiman_data/model/` for the files to be discoverable. New CaImAn users do not need to do anything as this is being taken care of during the installation process.
@@ -135,38 +144,9 @@ Alternative environments:
A list of known issues can be found [here](https://github.com/flatironinstitute/CaImAn/wiki/Known-Issues). If you still encounter problems please open an issue.
-
-## Graphical interface
-
-CaImAn comes with an experimental visual interface. To see an example on how use it, first load and run either
-* demo_OnACID_mesoscope.py
-* demo_caiman_basic.py
-* demo_pipeline.py
-
-Then you can start the visual interface by running the following command from the base caiman folder (make sure you are within your caiman environment):
-```
-ipython caiman/gui/gui_pyqtgraph_layout.py
-```
-You will then be asked to load the file that is generated at the end of each file (ending in .hdf5)
-
-A visual interface will appear, you will be able to:
-* regulate gain and contrast for the background image (correlation image)
-* regulate the threshold over the spatial masks to visualize components contours
-* click on neurons and see the corresponding trace and mask
-* select subset of neurons based on different quality metrics
-* save the resulting selection to a file in hdf5 format
-
-
-
-## Exporting results
-It is currently possible to export the results in two formats:
-* an HDF5 file that contains all the details of the CNMF object,
-the method to use is included in the
-* If you're using the [NWB](https://www.nwb.org/) file format, you can save the results of the analysis (`Estimates` object) back to the original NWB file. Consult this [demo](use_cases/NWB/demo_pipeline_NWB.py) for an example on how to use this feature.
-
## Documentation & Wiki
-Documentation of the code can be found [here](http://flatironinstitute.github.io/CaImAn/).
+Documentation of the code can be found [here](https://caiman.readthedocs.io/en/master/).
Moreover, our [wiki page](https://github.com/flatironinstitute/CaImAn/wiki) covers some aspects of the code.
### Installation for behavioral analysis
@@ -272,6 +252,7 @@ Special thanks to the following people for letting us use their datasets for our
* Weijian Yang, Darcy Peterka, Rafael Yuste, Columbia University
* Sue Ann Koay, David Tank, Princeton University
* Manolis Froudarakis, Jake Reimers, Andreas Tolias, Baylor College of Medicine
+* Clay Lacefield, Randy Bruno, Columbia University
## License
diff --git a/caiman/base/movies.py b/caiman/base/movies.py
index 56be4b447..4b47ef2c1 100644
--- a/caiman/base/movies.py
+++ b/caiman/base/movies.py
@@ -5,12 +5,6 @@
Contains the movie class.
-See Also:
-------------
-
-@url
-.. image::
-@author andrea giovannucci , deep-introspection
"""
# \package caiman/dource_ectraction/cnmf
@@ -83,7 +77,8 @@ class movie(ts.timeseries):
start_time=0
m=movie(input_arr, start_time=0,fr=33);
-
+ See https://docs.scipy.org/doc/numpy/user/basics.subclassing.html for
+ notes on objects that are descended from ndarray
"""
def __new__(cls, input_arr, **kwargs):
@@ -115,7 +110,7 @@ def apply_shifts_online(self, xy_shifts, save_base_name=None):
else:
return apply_shift_online(self, xy_shifts, save_base_name=save_base_name)
- def calc_min(self):
+ def calc_min(self) -> 'movie':
# todo: todocument
tmp = []
@@ -128,12 +123,11 @@ def calc_min(self):
return movie(input_arr=minval)
def motion_correct(self,
- max_shift_w=5,
- max_shift_h=5,
+ max_shift_w=5, max_shift_h=5,
num_frames_template=None,
template=None,
- method='opencv',
- remove_blanks=False, interpolation='cubic'):
+ method:str='opencv',
+ remove_blanks:bool=False, interpolation:str='cubic') -> Tuple[Any, Tuple, Any, Any]:
"""
Extract shifts and motion corrected movie automatically,
@@ -201,8 +195,8 @@ def motion_correct(self,
if remove_blanks:
max_h, max_w = np.max(shifts, axis=0)
min_h, min_w = np.min(shifts, axis=0)
- self = self.crop(crop_top=max_h, crop_bottom=-min_h + 1,
- crop_left=max_w, crop_right=-min_w, crop_begin=0, crop_end=0)
+ self.crop(crop_top=max_h, crop_bottom=-min_h + 1,
+ crop_left=max_w, crop_right=-min_w, crop_begin=0, crop_end=0)
return self, shifts, xcorrs, template
@@ -286,7 +280,8 @@ def motion_correct_3d(self,
return self, shifts, xcorrs, template
- def bin_median(self, window=10):
+
+ def bin_median(self, window:int=10) -> np.ndarray:
""" compute median of 3D array in along axis o by binning values
Args:
@@ -417,7 +412,7 @@ def extract_shifts(self, max_shift_w:int=5, max_shift_h:int=5, template=None, me
return (shifts, xcorrs)
- def apply_shifts(self, shifts, interpolation='linear', method='opencv', remove_blanks=False):
+ def apply_shifts(self, shifts, interpolation:str='linear', method:str='opencv', remove_blanks:bool=False):
"""
Apply precomputed shifts to a movie, using subpixels adjustment (cv2.INTER_CUBIC function)
@@ -426,6 +421,10 @@ def apply_shifts(self, shifts, interpolation='linear', method='opencv', remove_b
interpolation: 'linear', 'cubic', 'nearest' or cvs.INTER_XXX
+ method: (undocumented)
+
+ remove_blanks: (undocumented)
+
Returns:
self
@@ -498,8 +497,8 @@ def apply_shifts(self, shifts, interpolation='linear', method='opencv', remove_b
if remove_blanks:
max_h, max_w = np.max(shifts, axis=0)
min_h, min_w = np.min(shifts, axis=0)
- self = self.crop(crop_top=max_h, crop_bottom=-min_h + 1,
- crop_left=max_w, crop_right=-min_w, crop_begin=0, crop_end=0)
+ self.crop(crop_top=max_h, crop_bottom=-min_h + 1,
+ crop_left=max_w, crop_right=-min_w, crop_begin=0, crop_end=0)
return self
@@ -536,11 +535,17 @@ def linf(x, a, b):
return self
- def crop(self, crop_top=0, crop_bottom=0, crop_left=0, crop_right=0, crop_begin=0, crop_end=0):
- """ Crop movie
+ def crop(self, crop_top=0, crop_bottom=0, crop_left=0, crop_right=0, crop_begin=0, crop_end=0) -> None:
+ """
+ Crop movie (inline)
+
+ Args:
+ crop_top/crop_bottom/crop_left,crop_right: (undocumented)
+
+ crop_begin/crop_end: (undocumented)
"""
t, h, w = self.shape
- return self[crop_begin:t - crop_end, crop_top:h - crop_bottom, crop_left:w - crop_right]
+ self[:,:,:] = self[crop_begin:t - crop_end, crop_top:h - crop_bottom, crop_left:w - crop_right]
def computeDFF(self, secsWindow:int=5, quantilMin:int=8, method:str='only_baseline', order:str='F') -> Tuple[Any, Any]:
"""
@@ -1013,7 +1018,7 @@ def gaussian_blur_2D(self, kernel_size_x=5, kernel_size_y=5, kernel_std_x=1, ker
return self
- def median_blur_2D(self, kernel_size=3):
+ def median_blur_2D(self, kernel_size:float=3.0):
"""
Compute gaussian blut in 2D. Might be useful when motion correcting
@@ -1043,7 +1048,7 @@ def to_2D(self, order='F') -> np.ndarray:
d = d1 * d2
return np.reshape(self, (T, d), order=order)
- def zproject(self, method='mean', cmap=pl.cm.gray, aspect='auto', **kwargs) -> np.ndarray:
+ def zproject(self, method:str='mean', cmap=pl.cm.gray, aspect='auto', **kwargs) -> np.ndarray:
"""
Compute and plot projection across time:
@@ -1231,14 +1236,14 @@ def animate(i):
for i in range(10):
cv2.waitKey(100)
-def load(file_name, fr:float=30, start_time:float=0, meta_data:Dict=None, subindices=None,
+def load(file_name:Union[str,List[str]], fr:float=30, start_time:float=0, meta_data:Dict=None, subindices=None,
shape:Tuple[int,int]=None, var_name_hdf5:str='mov', in_memory:bool=False, is_behavior:bool=False,
bottom=0, top=0, left=0, right=0, channel=None, outtype=np.float32, is3D:bool=False) -> Any:
"""
load movie from file. Supports a variety of formats. tif, hdf5, npy and memory mapped. Matlab is experimental.
Args:
- file_name: string
+ file_name: string or List[str]
name of file. Possible extensions are tif, avi, npy, (npz and hdf5 are usable only if saved by calblitz)
fr: float
@@ -1285,7 +1290,7 @@ def load(file_name, fr:float=30, start_time:float=0, meta_data:Dict=None, subind
"""
# case we load movie from file
if max(top, bottom, left, right) > 0 and type(file_name) is str:
- file_name = [file_name]
+ file_name = [file_name] # type: ignore # mypy doesn't like that this changes type
if type(file_name) is list:
if shape is not None:
@@ -1346,7 +1351,7 @@ def load(file_name, fr:float=30, start_time:float=0, meta_data:Dict=None, subind
height = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT))
cv_failed = False
- dims:Union[Tuple,List] = [length, height, width]
+ dims = [length, height, width] # type: ignore # a list in one block and a tuple in another
if length == 0 or width == 0 or height == 0: #CV failed to load
cv_failed = True
if subindices is not None:
@@ -1485,7 +1490,7 @@ def rgb2gray(rgb):
elif extension == '.mmap':
filename = os.path.split(file_name)[-1]
- Yr, dims, T = load_memmap(os.path.join(
+ Yr, dims, T = load_memmap(os.path.join( # type: ignore # same dims typing issue as above
os.path.split(file_name)[0], filename))
images = np.reshape(Yr.T, [T] + list(dims), order='F')
if subindices is not None:
@@ -1652,7 +1657,7 @@ def sbxread(filename:str, k:int=0, n_frames=np.inf) -> np.ndarray:
# Paramters
N = max_idx + 1 # Last frame
- N = np.minimum(max_idx, n_frames)
+ N = np.minimum(N, n_frames)
nSamples = info['sz'][1] * info['recordsPerBuffer'] * 2 * info['nChan']
@@ -1665,8 +1670,11 @@ def sbxread(filename:str, k:int=0, n_frames=np.inf) -> np.ndarray:
x = ii16.max - np.fromfile(fo, dtype='uint16', count=int(nSamples / 2 * N))
x = x.reshape((int(info['nChan']), int(info['sz'][1]), int(
info['recordsPerBuffer']), int(N)), order='F')
+
x = x[0, :, :, :]
+ fo.close()
+
return x.transpose([2, 1, 0])
@@ -1682,7 +1690,6 @@ def sbxreadskip(filename:str, subindices:slice) -> np.ndarray:
if '.sbx' in filename:
filename = filename[:-4]
-
# Load info
info = loadmat_sbx(filename + '.mat')['info']
@@ -1707,18 +1714,17 @@ def sbxreadskip(filename:str, subindices:slice) -> np.ndarray:
start = 0
else:
start = subindices.start
-
+
if subindices.stop is None:
N = max_idx + 1 # Last frame
else:
- N = np.minimum(subindices.stop, max_idx + 1) .astype(np.int)
-
+ N = np.minimum(subindices.stop, max_idx + 1).astype(np.int)
+
if subindices.step is None:
skip = 1
else:
skip = subindices.step
-
iterable_elements = range(start, N, skip)
else:
@@ -1726,49 +1732,50 @@ def sbxreadskip(filename:str, subindices:slice) -> np.ndarray:
N = len(subindices)
iterable_elements = subindices
skip = 0
-
+
N_time = len(list(iterable_elements))
nSamples = info['sz'][1] * info['recordsPerBuffer'] * 2 * info['nChan']
+ assert nSamples >= 0
# Open File
fo = open(filename + '.sbx')
# Note: SBX files store the values strangely, its necessary to subtract the values from the max int16 to get the correct ones
-
-
-
+
counter = 0
-
+
if skip == 1:
# Note: SBX files store the values strangely, its necessary to subtract the values from the max int16 to get the correct ones
+ assert start * nSamples > 0
fo.seek(start * nSamples, 0)
ii16 = np.iinfo(np.uint16)
x = ii16.max - np.fromfile(fo, dtype='uint16', count=int(nSamples / 2 * (N-start)))
x = x.reshape((int(info['nChan']), int(info['sz'][1]), int(
info['recordsPerBuffer']), int(N-start)), order='F')
+
x = x[0, :, :, :]
-
+
else:
for k in iterable_elements:
- if counter%100 == 0:
- print('Reading Iteration:' + str(k))
+ assert k >= 0
+ if counter % 100 == 0:
+ print(f'Reading Iteration: {k}')
fo.seek(k * nSamples, 0)
ii16 = np.iinfo(np.uint16)
tmp = ii16.max - \
np.fromfile(fo, dtype='uint16', count=int(nSamples / 2 * 1))
-
+
tmp = tmp.reshape((int(info['nChan']), int(info['sz'][1]), int(
- info['recordsPerBuffer']), int(1)), order='F')
- if counter==0:
- x = np.zeros((tmp.shape[1],tmp.shape[2], tmp.shape[0],N_time))
- x[:,:,:,0] = tmp
- else:
- x[:,:,:,counter] = tmp
-
+ info['recordsPerBuffer'])), order='F')
+ if counter == 0:
+ x = np.zeros((tmp.shape[0], tmp.shape[1], tmp.shape[2], N_time))
+
+ x[:, :, :, counter] = tmp
counter += 1
- x = x[:, :, 0, :]
+ x = x[0, :, :, :]
+ fo.close()
return x.transpose([2, 1, 0])
diff --git a/caiman/base/rois.py b/caiman/base/rois.py
index e14866c51..0e6a2a7bb 100644
--- a/caiman/base/rois.py
+++ b/caiman/base/rois.py
@@ -32,7 +32,7 @@
from skimage.draw import polygon
import tempfile
import time
-from typing import List
+from typing import Any, Dict, List, Optional, Tuple
import zipfile
from ..motion_correction import tile_and_correct
@@ -42,7 +42,7 @@
except:
pass
-def com(A, d1, d2, d3=None):
+def com(A:np.ndarray, d1:int, d2:int, d3:Optional[int]=None) -> np.array:
"""Calculation of the center of mass for spatial components
Args:
@@ -78,7 +78,7 @@ def com(A, d1, d2, d3=None):
cm = (Coor * A / A.sum(axis=0)).T
return np.array(cm)
-def extract_binary_masks_from_structural_channel(Y, min_area_size=30, min_hole_size=15, gSig=5, expand_method='closing', selem=np.ones((3, 3))):
+def extract_binary_masks_from_structural_channel(Y, min_area_size:int=30, min_hole_size:int=15, gSig:int=5, expand_method:str='closing', selem:np.array=np.ones((3, 3))) -> Tuple[np.ndarray, np.array]:
"""Extract binary masks by using adaptive thresholding on a structural channel
Args:
@@ -142,7 +142,7 @@ def mask_to_2d(mask):
dims = np.shape(mask)
return scipy.sparse.coo_matrix(np.reshape(mask, (np.prod(dims), -1,), order='F'))
-def get_distance_from_A(masks_gt, masks_comp, min_dist=10):
+def get_distance_from_A(masks_gt, masks_comp, min_dist=10) -> List:
# todo todocument
_, d1, d2 = np.shape(masks_gt)
@@ -661,7 +661,7 @@ def norm_nrg(a_):
a[indx] = cumEn
return a.reshape(dims, order='F')
-def distance_masks(M_s, cm_s, max_dist, enclosed_thr=None):
+def distance_masks(M_s:List, cm_s:List[List], max_dist:float, enclosed_thr:Optional[float]=None) -> List:
"""
Compute distance matrix based on an intersection over union metric. Matrix are compared in order,
with matrix i compared with matrix i+1
@@ -745,7 +745,7 @@ def distance_masks(M_s, cm_s, max_dist, enclosed_thr=None):
D_s.append(D)
return D_s
-def find_matches(D_s, print_assignment=False):
+def find_matches(D_s, print_assignment:bool=False) -> Tuple[List, List]:
# todo todocument
matches = []
@@ -777,7 +777,7 @@ def find_matches(D_s, print_assignment=False):
# send back the results in the format we want
return matches, costs
-def link_neurons(matches, costs, max_cost=0.6, min_FOV_present=None):
+def link_neurons(matches:List[List[Tuple]], costs:List[List], max_cost:float=0.6, min_FOV_present:Optional[int]=None):
"""
Link neurons from different FOVs given matches and costs obtained from the hungarian algorithm
@@ -826,7 +826,7 @@ def link_neurons(matches, costs, max_cost=0.6, min_FOV_present=None):
logging.info(('num_neurons:' + str(num_neurons)))
return neurons
-def nf_load_masks(file_name, dims):
+def nf_load_masks(file_name:str, dims:Tuple[int, ...]) -> np.array:
# todo todocument
# load the regions (training data only)
@@ -841,7 +841,7 @@ def tomask(coords):
masks = np.array([tomask(s['coordinates']) for s in regions])
return masks
-def nf_masks_to_json(binary_masks, json_filename):
+def nf_masks_to_json(binary_masks:np.ndarray, json_filename:str) -> List[Dict]:
"""
Take as input a tensor of binary mask and produces json format for neurofinder
@@ -865,8 +865,7 @@ def nf_masks_to_json(binary_masks, json_filename):
return regions
-#%%
-def nf_masks_to_neurof_dict(binary_masks, dataset_name):
+def nf_masks_to_neurof_dict(binary_masks:np.ndarray, dataset_name:str) -> Dict[str, Any]:
"""
Take as input a tensor of binary mask and produces dict format for neurofinder
@@ -888,11 +887,7 @@ def nf_masks_to_neurof_dict(binary_masks, dataset_name):
return dset
-#%%
-
-
-
-def nf_read_roi(fileobj):
+def nf_read_roi(fileobj) -> np.ndarray:
'''
points = read_roi(fileobj)
Read ImageJ's ROI format
@@ -995,7 +990,7 @@ def getfloat():
return points
-def nf_read_roi_zip(fname, dims, return_names=False):
+def nf_read_roi_zip(fname:str, dims:Tuple[int, ...], return_names=False) -> np.array:
# todo todocument
with zipfile.ZipFile(fname) as zf:
@@ -1017,7 +1012,7 @@ def tomask(coords):
else:
return masks
-def nf_merge_roi_zip(fnames, idx_to_keep, new_fold):
+def nf_merge_roi_zip(fnames:List[str], idx_to_keep:List[List], new_fold:str):
"""
Create a zip file containing ROIs for ImageJ by combining elements from a list of ROI zip files
@@ -1054,13 +1049,13 @@ def nf_merge_roi_zip(fnames, idx_to_keep, new_fold):
shutil.make_archive(new_fold, 'zip', new_fold)
shutil.rmtree(new_fold)
-def extract_binary_masks_blob(A, neuron_radius, dims, num_std_threshold=1, minCircularity=0.5,
- minInertiaRatio=0.2, minConvexity=.8):
+def extract_binary_masks_blob(A, neuron_radius:float, dims:Tuple[int, ...], num_std_threshold:int=1, minCircularity:float=0.5,
+ minInertiaRatio:float=0.2, minConvexity:float=.8) -> Tuple[np.array, np.array, np.array]:
"""
Function to extract masks from data. It will also perform a preliminary selectino of good masks based on criteria like shape and size
Args:
- A: scipy.sparse matris
+ A: scipy.sparse matrix
contains the components as outputed from the CNMF algorithm
neuron_radius: float
@@ -1153,7 +1148,7 @@ def extract_binary_masks_blob(A, neuron_radius, dims, num_std_threshold=1, minCi
def extract_binary_masks_blob_parallel(A, neuron_radius, dims, num_std_threshold=1, minCircularity=0.5,
- minInertiaRatio=0.2, minConvexity=.8, dview=None):
+ minInertiaRatio=0.2, minConvexity=.8, dview=None) -> Tuple[List, List, List]:
# todo todocument
pars = []
@@ -1178,7 +1173,7 @@ def extract_binary_masks_blob_parallel(A, neuron_radius, dims, num_std_threshold
return masks, is_pos, is_neg
-def extract_binary_masks_blob_parallel_place_holder(pars):
+def extract_binary_masks_blob_parallel_place_holder(pars:Tuple) -> Tuple[Any, Any, Any]:
A, neuron_radius, dims, num_std_threshold, _, minInertiaRatio, minConvexity = pars
masks_ws, pos_examples, neg_examples = extract_binary_masks_blob(A, neuron_radius,
dims, num_std_threshold=num_std_threshold, minCircularity=0.5,
@@ -1187,14 +1182,14 @@ def extract_binary_masks_blob_parallel_place_holder(pars):
-def extractROIsFromPCAICA(spcomps, numSTD=4, gaussiansigmax=2, gaussiansigmay=2, thresh=None):
+def extractROIsFromPCAICA(spcomps, numSTD=4, gaussiansigmax=2, gaussiansigmay=2, thresh=None) -> Tuple[List[np.array], List]:
"""
Given the spatial components output of the IPCA_stICA function extract possible regions of interest
The algorithm estimates the significance of a components by thresholding the components after gaussian smoothing
Args:
- spcompomps, 3d array containing the spatial components
+ spcomps: 3d array containing the spatial components
numSTD: number of standard deviation above the mean of the spatial component to be considered signiificant
"""
@@ -1227,8 +1222,8 @@ def extractROIsFromPCAICA(spcomps, numSTD=4, gaussiansigmax=2, gaussiansigmay=2,
return allMasks, maskgrouped
-def detect_duplicates_and_subsets(binary_masks, predictions=None, r_values=None, dist_thr=0.1, min_dist=10,
- thresh_subset=0.8):
+def detect_duplicates_and_subsets(binary_masks, predictions=None, r_values=None, dist_thr:float=0.1, min_dist=10,
+ thresh_subset:float=0.8):
cm = [scipy.ndimage.center_of_mass(mm) for mm in binary_masks]
sp_rois = scipy.sparse.csc_matrix(
@@ -1319,7 +1314,7 @@ def detect_duplicates_and_subsets(binary_masks, predictions=None, r_values=None,
return indices_orig, indices_to_keep, indices_to_remove, D, overlap
-def detect_duplicates(file_name, dist_thr=0.1, FOV=(512, 512)):
+def detect_duplicates(file_name:str, dist_thr:float=0.1, FOV:Tuple[int, ...]=(512, 512)) -> Tuple[List, List]:
"""
Removes duplicate ROIs from file file_name
diff --git a/caiman/base/timeseries.py b/caiman/base/timeseries.py
index fe805c44d..29b09c8d3 100644
--- a/caiman/base/timeseries.py
+++ b/caiman/base/timeseries.py
@@ -129,15 +129,14 @@ def save(self,
var_name_hdf5='mov',
sess_desc='some_description',
identifier='some identifier',
- exp_desc='experiment description',
imaging_plane_description='some imaging plane description',
emission_lambda=520.0,
indicator='OGB-1',
location='brain',
starting_time=0.,
experimenter='Dr Who',
- lab_name='',
- institution='',
+ lab_name=None,
+ institution=None,
experiment_description='Experiment Description',
session_id='Session ID'):
"""
diff --git a/caiman/behavior/behavior.py b/caiman/behavior/behavior.py
index 889031c33..254ad7ce4 100644
--- a/caiman/behavior/behavior.py
+++ b/caiman/behavior/behavior.py
@@ -9,27 +9,28 @@
@author: agiovann
"""
-#%%
from builtins import zip
from builtins import range
-import time
from past.utils import old_div
-import caiman as cm
+
+import cv2
import numpy as np
import pylab as pl
import scipy
from scipy.sparse import coo_matrix
from scipy.io import loadmat
-import cv2
from sklearn.decomposition import NMF
+import time
+from typing import List
+
+import caiman as cm
try:
cv2.setNumThreads(0)
except:
pass
-#%% dense flow
-def select_roi(img, n_rois=1):
+def select_roi(img:np.ndarray, n_rois:int=1) -> List:
"""
Create a mask from a the convex polygon enclosed between selected points
@@ -58,23 +59,16 @@ def select_roi(img, n_rois=1):
return masks
-#%%
def to_polar(x, y):
mag, ang = cv2.cartToPolar(x, y)
return mag, ang
-#%%
-
-
def get_nonzero_subarray(arr, mask):
x, y = mask.nonzero()
return arr.toarray()[x.min():x.max() + 1, y.min():y.max() + 1]
-#%%
-
-
-def extract_motor_components_OF(m, n_components, mask=None, resize_fact=.5, only_magnitude=False, max_iter=1000,
- verbose=False, method_factorization='nmf', max_iter_DL=-30):
+def extract_motor_components_OF(m, n_components, mask=None, resize_fact:float=.5, only_magnitude:bool=False, max_iter:int=1000,
+ verbose:bool=False, method_factorization:str='nmf', max_iter_DL=-30) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
# todo todocument
if mask is not None:
mask = coo_matrix(np.array(mask).squeeze())
@@ -104,10 +98,8 @@ def extract_motor_components_OF(m, n_components, mask=None, resize_fact=.5, only
return spatial_filter_, time_trace_, of_or
-
-#%%
def extract_magnitude_and_angle_from_OF(spatial_filter_, time_trace_, of_or,
- num_std_mag_for_angle=.6, sav_filter_size=3, only_magnitude=False):
+ num_std_mag_for_angle=.6, sav_filter_size=3, only_magnitude=False) -> Tuple[List, List, List, List]:
# todo todocument
mags = []
@@ -150,10 +142,8 @@ def extract_magnitude_and_angle_from_OF(spatial_filter_, time_trace_, of_or,
return mags, dircts, dircts_thresh, spatial_masks_thr
-
-#%%
-def compute_optical_flow(m, mask=None, polar_coord=True, do_show=False, do_write=False, file_name=None, gain_of=None,
- frate=30, pyr_scale=.1, levels=3, winsize=25, iterations=3, poly_n=7, poly_sigma=1.5):
+def compute_optical_flow(m, mask=None, polar_coord:bool=True, do_show:bool=False, do_write:bool=False, file_name=None, gain_of=None,
+ frate:float=30.0, pyr_scale:float=.1, levels:int=3, winsize:int=25, iterations:int=3, poly_n=7, poly_sigma=1.5):
"""
This function compute the optical flow of behavioral movies using the opencv cv2.calcOpticalFlowFarneback function
@@ -264,7 +254,7 @@ def compute_optical_flow(m, mask=None, polar_coord=True, do_show=False, do_write
#%% NMF
-def extract_components(mov_tot, n_components=6, normalize_std=True, max_iter_DL=-30, method_factorization='nmf', **kwargs):
+def extract_components(mov_tot, n_components:int=6, normalize_std:bool=True, max_iter_DL=-30, method_factorization:str='nmf', **kwargs) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
From optical flow images can extract spatial and temporal components
@@ -286,7 +276,7 @@ def extract_components(mov_tot, n_components=6, normalize_std=True, max_iter_DL=
spatial_filter: ndarray
set of spatial inferred filters
- time_trace:ndarray
+ time_trace: ndarray
set of time components
norm_fact: ndarray
@@ -336,7 +326,7 @@ def extract_components(mov_tot, n_components=6, normalize_std=True, max_iter_DL=
print(el_t)
return spatial_filter, time_trace, norm_fact
-def plot_components(sp_filt, t_trace):
+def plot_components(sp_filt, t_trace) -> None:
# todo: todocument
pl.figure()
count = 0
diff --git a/caiman/cluster.py b/caiman/cluster.py
index ba83c5fa6..2eb97a39c 100644
--- a/caiman/cluster.py
+++ b/caiman/cluster.py
@@ -1,6 +1,5 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
-
""" functions related to the creation and management of the cluster
We put arrays on disk as raw bytes, extending along the first dimension.
@@ -40,7 +39,12 @@
logger = logging.getLogger(__name__)
-def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[List[int],Tuple], border_pix:int=0, indices=[slice(None)]*2) -> Tuple[List, List]:
+
+def extract_patch_coordinates(dims: Tuple,
+ rf: Union[List, Tuple],
+ stride: Union[List[int], Tuple],
+ border_pix: int = 0,
+ indices=[slice(None)] * 2) -> Tuple[List, List]:
"""
Partition the FOV in patches
and return the indexed in 2D and 1D (flatten, order='F') formats
@@ -58,25 +62,24 @@ def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[Lis
sl_start = [0 if sl.start is None else sl.start for sl in indices]
sl_stop = [dim if sl.stop is None else sl.stop for (sl, dim) in zip(indices, dims)]
- sl_step = [1 for sl in indices] # not used
+ sl_step = [1 for sl in indices] # not used
dims_large = dims
dims = np.minimum(np.array(dims) - border_pix, sl_stop) - np.maximum(border_pix, sl_start)
coords_flat = []
shapes = []
- iters = [list(range(rf[i], dims[i] - rf[i], 2 * rf[i] - stride[i])) + [dims[i] - rf[i]]
- for i in range(len(dims))]
+ iters = [list(range(rf[i], dims[i] - rf[i], 2 * rf[i] - stride[i])) + [dims[i] - rf[i]] for i in range(len(dims))]
coords = np.empty(list(map(len, iters)) + [len(dims)], dtype=np.object)
for count_0, xx in enumerate(iters[0]):
coords_x = np.arange(xx - rf[0], xx + rf[0] + 1)
coords_x = coords_x[(coords_x >= 0) & (coords_x < dims[0])]
- coords_x += border_pix*0 + np.maximum(sl_start[0], border_pix)
+ coords_x += border_pix * 0 + np.maximum(sl_start[0], border_pix)
for count_1, yy in enumerate(iters[1]):
coords_y = np.arange(yy - rf[1], yy + rf[1] + 1)
coords_y = coords_y[(coords_y >= 0) & (coords_y < dims[1])]
- coords_y += border_pix*0 + np.maximum(sl_start[1], border_pix)
+ coords_y += border_pix * 0 + np.maximum(sl_start[1], border_pix)
if len(dims) == 2:
idxs = np.meshgrid(coords_x, coords_y)
@@ -86,7 +89,7 @@ def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[Lis
coords_ = np.ravel_multi_index(idxs, dims_large, order='F')
coords_flat.append(coords_.flatten())
- else: # 3D data
+ else: # 3D data
if border_pix > 0:
raise Exception(
@@ -109,7 +112,8 @@ def extract_patch_coordinates(dims:Tuple, rf:Union[List,Tuple], stride:Union[Lis
#%%
-def apply_to_patch(mmap_file, shape:Tuple[Any,Any,Any], dview, rf, stride, function, *args, **kwargs) -> Tuple[List,Any,Tuple]:
+def apply_to_patch(mmap_file, shape: Tuple[Any, Any, Any], dview, rf, stride, function, *args,
+ **kwargs) -> Tuple[List, Any, Tuple]:
"""
apply function to patches in parallel or not
@@ -152,11 +156,9 @@ def apply_to_patch(mmap_file, shape:Tuple[Any,Any,Any], dview, rf, stride, funct
stride1 = stride
stride2 = stride
- idx_flat, idx_2d = extract_patch_coordinates(
- (d1, d2), rf=(rf1, rf2), stride=(stride1, stride2))
+ idx_flat, idx_2d = extract_patch_coordinates((d1, d2), rf=(rf1, rf2), stride=(stride1, stride2))
- shape_grid = tuple(np.ceil(
- (d1 * 1. / (rf1 * 2 - stride1), d2 * 1. / (rf2 * 2 - stride2))).astype(np.int))
+ shape_grid = tuple(np.ceil((d1 * 1. / (rf1 * 2 - stride1), d2 * 1. / (rf2 * 2 - stride2))).astype(np.int))
if d1 <= rf1 * 2:
shape_grid = (1, shape_grid[1])
if d2 <= rf2 * 2:
@@ -167,8 +169,7 @@ def apply_to_patch(mmap_file, shape:Tuple[Any,Any,Any], dview, rf, stride, funct
args_in = []
for id_f, id_2d in zip(idx_flat[:], idx_2d[:]):
- args_in.append((mmap_file.filename, id_f,
- id_2d, function, args, kwargs))
+ args_in.append((mmap_file.filename, id_f, id_2d, function, args, kwargs))
logger.debug("Flat index is of length " + str(len(idx_flat)))
if dview is not None:
@@ -184,10 +185,12 @@ def apply_to_patch(mmap_file, shape:Tuple[Any,Any,Any], dview, rf, stride, funct
file_res = list(map(function_place_holder, args_in))
return file_res, idx_flat, shape_grid
+
+
#%%
-def function_place_holder(args_in):
+def function_place_holder(args_in: Tuple) -> np.ndarray:
#todo: todocument
file_name, idx_, shapes, function, args, kwargs = args_in
@@ -195,8 +198,7 @@ def function_place_holder(args_in):
Yr = Yr[idx_, :]
Yr.filename = file_name
_, T = Yr.shape
- Y = np.reshape(Yr, (shapes[1], shapes[0], T),
- order='F').transpose([2, 0, 1])
+ Y = np.reshape(Yr, (shapes[1], shapes[0], T), order='F').transpose([2, 0, 1])
[T, d1, d2] = Y.shape
res_fun = function(Y, *args, **kwargs)
@@ -208,10 +210,11 @@ def function_place_holder(args_in):
return res_fun
+
#%%
-def start_server(slurm_script:str=None, ipcluster:str="ipcluster", ncpus:int=None) -> None:
+def start_server(slurm_script: str = None, ipcluster: str = "ipcluster", ncpus: int = None) -> None:
"""
programmatically start the ipyparallel server
@@ -220,7 +223,7 @@ def start_server(slurm_script:str=None, ipcluster:str="ipcluster", ncpus:int=Non
number of processors
ipcluster : str
- ipcluster binary file name; requires 4 path separators on Windows. ipcluster="C:\\\\Anaconda2\\\\Scripts\\\\ipcluster.exe"
+ ipcluster binary file name; requires 4 path separators on Windows. ipcluster="C:\\\\Anaconda3\\\\Scripts\\\\ipcluster.exe"
Default: "ipcluster"
"""
logger.info("Starting cluster...")
@@ -230,21 +233,21 @@ def start_server(slurm_script:str=None, ipcluster:str="ipcluster", ncpus:int=Non
if slurm_script is None:
if ipcluster == "ipcluster":
- subprocess.Popen(
- "ipcluster start -n {0}".format(ncpus), shell=True, close_fds=(os.name != 'nt'))
+ subprocess.Popen("ipcluster start -n {0}".format(ncpus), shell=True, close_fds=(os.name != 'nt'))
else:
- subprocess.Popen(shlex.split(
- "{0} start -n {1}".format(ipcluster, ncpus)), shell=True, close_fds=(os.name != 'nt'))
+ subprocess.Popen(shlex.split("{0} start -n {1}".format(ipcluster, ncpus)),
+ shell=True,
+ close_fds=(os.name != 'nt'))
time.sleep(1.5)
# Check that all processes have started
client = ipyparallel.Client()
time.sleep(1.5)
while len(client) < ncpus:
- sys.stdout.write(".") # Give some visual feedback of things starting
- sys.stdout.flush() # (de-buffered)
+ sys.stdout.write(".") # Give some visual feedback of things starting
+ sys.stdout.flush() # (de-buffered)
time.sleep(0.5)
logger.debug('Making sure everything is up and running')
- client.direct_view().execute('__a=1', block=True) # when done on all, we're set to go
+ client.direct_view().execute('__a=1', block=True) # when done on all, we're set to go
else:
shell_source(slurm_script)
pdir, profile = os.environ['IPPPDIR'], os.environ['IPPPROFILE']
@@ -256,18 +259,18 @@ def start_server(slurm_script:str=None, ipcluster:str="ipcluster", ncpus:int=Non
c.close()
sys.stdout.write("start_server: done\n")
-def shell_source(script:str) -> None:
+
+def shell_source(script: str) -> None:
""" Run a source-style bash script, copy resulting env vars to current process. """
# XXX This function is weird and maybe not a good idea. People easily might expect
# it to handle conditionals. Maybe just make them provide a key-value file
#introduce echo to indicate the end of the output
- pipe = subprocess.Popen(". %s; env; echo 'FINISHED_CLUSTER'" %
- script, stdout=subprocess.PIPE, shell=True)
+ pipe = subprocess.Popen(". %s; env; echo 'FINISHED_CLUSTER'" % script, stdout=subprocess.PIPE, shell=True)
env = dict()
while True:
line = pipe.stdout.readline().decode('utf-8').rstrip()
- if 'FINISHED_CLUSTER' in line: # find the keyword set above to determine the end of the output stream
+ if 'FINISHED_CLUSTER' in line: # find the keyword set above to determine the end of the output stream
break
logger.debug("shell_source parsing line[" + str(line) + "]")
lsp = str(line).split("=", 1)
@@ -277,7 +280,8 @@ def shell_source(script:str) -> None:
os.environ.update(env)
pipe.stdout.close()
-def stop_server(ipcluster:str='ipcluster', pdir:str=None, profile:str=None, dview=None) -> None:
+
+def stop_server(ipcluster: str = 'ipcluster', pdir: str = None, profile: str = None, dview=None) -> None:
"""
programmatically stops the ipyparallel server
@@ -315,7 +319,7 @@ def stop_server(ipcluster:str='ipcluster', pdir:str=None, profile:str=None, dvie
try:
shutil.rmtree('./log/')
except:
- logger.info('creating log folder') # FIXME Not what this means
+ logger.info('creating log folder') # FIXME Not what this means
files = glob.glob('*.log')
os.mkdir('./log')
@@ -325,11 +329,15 @@ def stop_server(ipcluster:str='ipcluster', pdir:str=None, profile:str=None, dvie
else:
if ipcluster == "ipcluster":
- proc = subprocess.Popen(
- "ipcluster stop", shell=True, stderr=subprocess.PIPE, close_fds=(os.name != 'nt'))
+ proc = subprocess.Popen("ipcluster stop",
+ shell=True,
+ stderr=subprocess.PIPE,
+ close_fds=(os.name != 'nt'))
else:
proc = subprocess.Popen(shlex.split(ipcluster + " stop"),
- shell=True, stderr=subprocess.PIPE, close_fds=(os.name != 'nt'))
+ shell=True,
+ stderr=subprocess.PIPE,
+ close_fds=(os.name != 'nt'))
line_out = proc.stderr.readline()
if b'CRITICAL' in line_out:
@@ -348,10 +356,15 @@ def stop_server(ipcluster:str='ipcluster', pdir:str=None, profile:str=None, dvie
proc.stderr.close()
logger.info("stop_cluster(): done")
+
+
#%%
-def setup_cluster(backend:str='multiprocessing', n_processes:int=None, single_thread:bool=False, ignore_preexisting:bool=False) -> Tuple[Any, Any, Optional[int]]:
+def setup_cluster(backend: str = 'multiprocessing',
+ n_processes: int = None,
+ single_thread: bool = False,
+ ignore_preexisting: bool = False) -> Tuple[Any, Any, Optional[int]]:
"""Setup and/or restart a parallel cluster.
Args:
backend: str
@@ -386,7 +399,7 @@ def setup_cluster(backend:str='multiprocessing', n_processes:int=None, single_th
stop_server()
except:
logger.debug('Nothing to stop')
- slurm_script = '/mnt/home/agiovann/SOFTWARE/CaImAn/SLURM/slurmStart.sh'
+ slurm_script = '/mnt/home/agiovann/SOFTWARE/CaImAn/SLURM/slurmStart.sh' # FIXME: Make this a documented environment variable
logger.info([str(n_processes), slurm_script])
start_server(slurm_script=slurm_script, ncpus=n_processes)
pdir, profile = os.environ['IPPPDIR'], os.environ['IPPPROFILE']
@@ -397,7 +410,7 @@ def setup_cluster(backend:str='multiprocessing', n_processes:int=None, single_th
stop_server()
start_server(ncpus=n_processes)
c = Client()
- logger.info('Started ipyparallel cluster: Using ' + str(len(c)) + ' processes')
+ logger.info(f'Started ipyparallel cluster: Using {len(c)} processes')
dview = c[:len(c)]
elif (backend == 'multiprocessing') or (backend == 'local'):
@@ -411,16 +424,16 @@ def setup_cluster(backend:str='multiprocessing', n_processes:int=None, single_th
'A cluster is already runnning. Terminate with dview.terminate() if you want to restart.')
if (platform.system() == 'Darwin') and (sys.version_info > (3, 0)):
try:
- if 'kernel' in get_ipython().trait_names(): # type: ignore
- # If you're on OSX and you're running under Jupyter or Spyder,
- # which already run the code in a forkserver-friendly way, this
- # can eliminate some setup and make this a reasonable approach.
- # Otherwise, seting VECLIB_MAXIMUM_THREADS=1 or using a different
- # blas/lapack is the way to avoid the issues.
- # See https://github.com/flatironinstitute/CaImAn/issues/206 for more
- # info on why we're doing this (for now).
+ if 'kernel' in get_ipython().trait_names(): # type: ignore
+ # If you're on OSX and you're running under Jupyter or Spyder,
+ # which already run the code in a forkserver-friendly way, this
+ # can eliminate some setup and make this a reasonable approach.
+ # Otherwise, seting VECLIB_MAXIMUM_THREADS=1 or using a different
+ # blas/lapack is the way to avoid the issues.
+ # See https://github.com/flatironinstitute/CaImAn/issues/206 for more
+ # info on why we're doing this (for now).
multiprocessing.set_start_method('forkserver', force=True)
- except: # If we're not running under ipython, don't do anything.
+ except: # If we're not running under ipython, don't do anything.
pass
c = None
diff --git a/caiman/components_evaluation.py b/caiman/components_evaluation.py
index 3d60dcdd9..30896ad7c 100644
--- a/caiman/components_evaluation.py
+++ b/caiman/components_evaluation.py
@@ -15,6 +15,7 @@
import logging
import numpy as np
import os
+import tensorflow as tf
import scipy
from scipy.sparse import csc_matrix
from scipy.stats import norm
@@ -211,7 +212,7 @@ def classify_components_ep(Y, A, C, b, f, Athresh=0.1, Npeaks=5, tB=-3, tA=10, t
significant_samples:List[Any] = []
for i in range(K):
- if i % 200 == 0: # Show status periodically
+ if (i+1) % 200 == 0: # Show status periodically
logging.info('Components evaluated:' + str(i))
if LOC[i] is not None:
atemp = A[:, i].toarray().flatten()
@@ -257,17 +258,15 @@ def evaluate_components_CNN(A, dims, gSig, model_name:str=os.path.join(caiman_da
import os
if not isGPU:
-
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
try:
os.environ["KERAS_BACKEND"] = "tensorflow"
- from keras.models import model_from_json
+ from tensorflow.keras.models import model_from_json
use_keras = True
- logging.debug('Using Keras')
+ logging.info('Using Keras')
except(ModuleNotFoundError):
- import tensorflow as tf
use_keras = False
- logging.debug('Using Tensorflow')
+ logging.info('Using Tensorflow')
if loaded_model is None:
if use_keras:
@@ -295,7 +294,8 @@ def evaluate_components_CNN(A, dims, gSig, model_name:str=os.path.join(caiman_da
raise FileNotFoundError("File for requested model {} not found".format(model_name))
loaded_model = load_graph(model_file)
- logging.info("Loaded model from disk")
+ logging.debug("Loaded model from disk")
+
half_crop = np.minimum(
gSig[0] * 4 + 1, patch_size), np.minimum(gSig[1] * 4 + 1, patch_size)
dims = np.array(dims)
@@ -317,6 +317,7 @@ def evaluate_components_CNN(A, dims, gSig, model_name:str=os.path.join(caiman_da
with tf.Session(graph=loaded_model) as sess:
predictions = sess.run(
tf_out, feed_dict={tf_in: final_crops[:, :, :, np.newaxis]})
+ sess.close()
return predictions, final_crops
#%%
@@ -326,11 +327,11 @@ def evaluate_components(Y:np.ndarray, traces:np.ndarray, A, C, b, f, final_frate
Athresh:float=0.1, Npeaks:int=5, thresh_C:float=0.3, sigma_factor:float=3.) -> Tuple[Any, Any, Any, Any, Any, Any]:
""" Define a metric and order components according to the probability of some "exceptional events" (like a spike).
- Such probability is defined as the likeihood of observing the actual trace value over N samples given an estimated noise distribution.
+ Such probability is defined as the likelihood of observing the actual trace value over N samples given an estimated noise distribution.
The function first estimates the noise distribution by considering the dispersion around the mode.
This is done only using values lower than the mode.
The estimation of the noise std is made robust by using the approximation std=iqr/1.349.
- Then, the probavility of having N consecutive eventsis estimated.
+ Then, the probability of having N consecutive events is estimated.
This probability is used to order the components.
The algorithm also measures the reliability of the spatial mask by comparing the filters in A
with the average of the movies over samples where exceptional events happen, after removing (if possible)
@@ -397,7 +398,7 @@ def evaluate_components(Y:np.ndarray, traces:np.ndarray, A, C, b, f, final_frate
Yr = np.reshape(Y, (np.prod(dims), T), order='F')
- logging.info('Computing event exceptionality delta')
+ logging.debug('Computing event exceptionality delta')
fitness_delta, erfc_delta, _, _ = compute_event_exceptionality(
np.diff(traces, axis=1), robust_std=robust_std, N=N, sigma_factor=sigma_factor)
@@ -435,11 +436,11 @@ def evaluate_components(Y:np.ndarray, traces:np.ndarray, A, C, b, f, final_frate
else:
traces -= tr_BL[padbefore:-padafter].T
- logging.info('Computing event exceptionality')
+ logging.debug('Computing event exceptionality')
fitness_raw, erfc_raw, _, _ = compute_event_exceptionality(
traces, robust_std=robust_std, N=N, sigma_factor=sigma_factor)
- logging.info('Evaluating spatial footprint')
+ logging.debug('Evaluating spatial footprint')
# compute the overlap between spatial and movie average across samples with significant events
r_values, significant_samples = classify_components_ep(Yr, A, C, b, f, Athresh=Athresh, Npeaks=Npeaks, tB=tB,
tA=tA, thres=thresh_C)
@@ -539,7 +540,6 @@ def estimate_components_quality_auto(Y, A, C, b, f, YrA, frate, decay_time, gSig
-min_SNR_reject) * N_samples
traces = C + YrA
-
_, _, fitness_raw, _, r_values = estimate_components_quality( # type: ignore # mypy cannot reason about return_all
traces, Y, A, C, b, f, final_frate=frate, Npeaks=Npeaks, r_values_min=r_values_min, fitness_min=fitness_min,
fitness_delta_min=thresh_fitness_delta, return_all=True, dview=dview, num_traces_per_group=50, N=N_samples)
@@ -610,11 +610,11 @@ def estimate_components_quality(traces, Y, A, C, b, f, final_frate=30, Npeaks=10
remove_baseline=True, dview=None, robust_std=False, Athresh=0.1, thresh_C=0.3, num_traces_per_group=20) -> Tuple[np.ndarray, ...]:
""" Define a metric and order components according to the probability of some "exceptional events" (like a spike).
- Such probability is defined as the likeihood of observing the actual trace value over N samples given an estimated noise distribution.
+ Such probability is defined as the likelihood of observing the actual trace value over N samples given an estimated noise distribution.
The function first estimates the noise distribution by considering the dispersion around the mode.
This is done only using values lower than the mode.
The estimation of the noise std is made robust by using the approximation std=iqr/1.349.
- Then, the probavility of having N consecutive eventsis estimated.
+ Then, the probability of having N consecutive events is estimated.
This probability is used to order the components.
The algorithm also measures the reliability of the spatial mask by comparing the filters in A
with the average of the movies over samples where exceptional events happen, after removing (if possible)
@@ -693,7 +693,7 @@ def estimate_components_quality(traces, Y, A, C, b, f, final_frate=30, Npeaks=10
if dview is None:
res = map(evaluate_components_placeholder, params)
else:
- logging.info('EVALUATING IN PARALLEL... NOT RETURNING ERFCs')
+ logging.info('Component evaluation in parallel')
if 'multiprocessing' in str(type(dview)):
res = dview.map_async(
evaluate_components_placeholder, params).get(4294967)
diff --git a/caiman/gui/gui_pyqtgraph_layout.py b/caiman/gui/gui_pyqtgraph_layout.py
index 0ec1d5518..53030c3a8 100644
--- a/caiman/gui/gui_pyqtgraph_layout.py
+++ b/caiman/gui/gui_pyqtgraph_layout.py
@@ -1,15 +1,15 @@
+import cv2
+import numpy as np
import pyqtgraph as pg
+import scipy
from pyqtgraph import FileDialog
-from pyqtgraph.Qt import QtCore, QtGui
-from pyqtgraph.parametertree import Parameter, ParameterTree, ParameterItem, registerParameterType
-import numpy as np
+from pyqtgraph.Qt import QtGui
+from pyqtgraph.parametertree import Parameter, ParameterTree
+
import caiman as cm
from caiman.source_extraction.cnmf.cnmf import load_CNMF
-from caiman.source_extraction.cnmf.online_cnmf import load_OnlineCNMF
+from caiman.source_extraction.cnmf.params import CNMFParams
-import cv2
-import scipy
-import os
# Always start by initializing Qt (only once per application)
app = QtGui.QApplication([])
@@ -28,74 +28,112 @@
except NameError:
print('Not launched under iPython')
-def make_color_img(img, gain=255, out_type=np.uint8):
- min = img.min()
- max = img.max()
- img = (img-min)/(max-min)*gain
+
+def make_color_img(img, gain=255, min_max=None,out_type=np.uint8):
+ if min_max is None:
+ min_ = img.min()
+ max_ = img.max()
+ else:
+ min_, max_ = min_max
+
+ img = (img-min_)/(max_-min_)*gain
img = img.astype(out_type)
img = np.dstack([img]*3)
return img
+
F = FileDialog()
# load object saved by CNMF
# cnm_obj = load_CNMF('/Users/agiovann/caiman_data/example_movies/memmap__d1_60_d2_80_d3_1_order_C_frames_2000_save.hdf5')
-cnm_obj = load_CNMF(F.getOpenFileName(caption='Load CNMF Object',filter='*.hdf5')[0])
-
+fpath = F.getOpenFileName(caption='Load CNMF Object', filter='HDF5 (*.h5 *.hdf5);;NWB (*.nwb)')[0]
+cnm_obj = load_CNMF(fpath)
# movie
# mov = cm.load('/Users/agiovann/caiman_data/example_movies/memmap__d1_60_d2_80_d3_1_order_C_frames_2000_.mmap')
-mov = cm.load(cnm_obj.mmap_file)
-
+if fpath[-3:] == 'nwb':
+ mov = cm.load(cnm_obj.mmap_file, var_name_hdf5='acquisition/TwoPhotonSeries')
+else:
+ mov = cm.load(cnm_obj.mmap_file)
# load summary image
# Cn = cm.load('/Users/agiovann/caiman_data/example_movies/memmap__d1_60_d2_80_d3_1_order_C_frames_2000__Cn.tif')
-Cn = cnm_obj.estimates.Cn
+estimates = cnm_obj.estimates
+params_obj = cnm_obj.params
+
+if not hasattr(estimates, 'Cn'):
+ estimates.Cn = cm.local_correlations(mov, swap_dim=False)
+
+min_mov = np.min(mov)
+max_mov = np.max(mov)
+
+min_mov_denoise = np.min(estimates.A.dot(estimates.C))
+max_mov_denoise = np.max(estimates.A.dot(estimates.C))
+background_num = -1
+neuron_selected = False
+nr_index = 0
+
+min_background = np.min(estimates.b, axis=0)*np.min(estimates.f, axis=1)
+max_background = np.max(estimates.b, axis=0)*np.max(estimates.f, axis=1)
-estimates = cnm_obj.estimates
if not hasattr(estimates, 'accepted_list'):
# if estimates.discarded_components.A.shape[-1] > 0:
# estimates.restore_discarded_components()
estimates.accepted_list = np.array([], dtype=np.int)
estimates.rejected_list = np.array([], dtype=np.int)
- estimates.img_components = estimates.A.toarray().reshape((estimates.dims[0], estimates.dims[1],-1), order='F').transpose([2,0,1])
+ estimates.img_components = estimates.A.toarray().reshape((estimates.dims[0], estimates.dims[1], -1), order='F').transpose([2,0,1])
estimates.cms = np.array([scipy.ndimage.measurements.center_of_mass(comp) for comp in estimates.img_components])
estimates.idx_components = np.arange(estimates.nr)
estimates.idx_components_bad = np.array([])
- estimates.background_image = make_color_img(Cn)
+ estimates.background_image = make_color_img(estimates.Cn)
# Generate image data
- estimates.img_components /= estimates.img_components.max(axis=(1,2))[:,None,None]
+ estimates.img_components /= estimates.img_components.max(axis=(1, 2))[:, None, None]
estimates.img_components *= 255
estimates.img_components = estimates.img_components.astype(np.uint8)
+def draw_contours_overall(md):
+ if md is "reset":
+ draw_contours()
+ elif md is "neurons":
+ if neuron_selected is True:
+ #if a specific neuron has been selected, only one contour should be changed while thrshcomp_line is changing
+ if nr_index is 0:
+ #if user does not start to move through the frames
+ draw_contours_update(estimates.background_image, img)
+ draw_contours_update(comp2_scaled, img2)
+ else:
+ draw_contours_update(raw_mov_scaled, img)
+ draw_contours_update(frame_denoise_scaled, img2)
+ else:
+ #if no specific neuron has been selected, all the contours are changing
+ draw_contours()
+ else:
+ #md is "background":
+ return
+
-
-#%%
def draw_contours():
- global thrshcomp_line, estimates, cnm_obj, img
+ global thrshcomp_line, estimates, img
bkgr_contours = estimates.background_image.copy()
if len(estimates.idx_components) > 0:
contours = [cv2.findContours(cv2.threshold(img, np.int(thrshcomp_line.value()), 255, 0)[1], cv2.RETR_TREE,
cv2.CHAIN_APPROX_SIMPLE)[0] for img in estimates.img_components[estimates.idx_components]]
-
SNRs = np.array(estimates.r_values)
iidd = np.array(estimates.idx_components)
- idx1 = np.where(SNRs[iidd]<0.1)[0]
- idx2 = np.where((SNRs[iidd]>=0.1) &
- (SNRs[iidd]<0.25))[0]
- idx3 = np.where((SNRs[iidd]>=0.25) &
- (SNRs[iidd]<0.5))[0]
- idx4 = np.where((SNRs[iidd]>=0.5) &
- (SNRs[iidd]<0.75))[0]
- idx5 = np.where((SNRs[iidd]>=0.75) &
- (SNRs[iidd]<0.9))[0]
- idx6 = np.where(SNRs[iidd]>=0.9)[0]
-
-
-
+ idx1 = np.where(SNRs[iidd] < 0.1)[0]
+ idx2 = np.where((SNRs[iidd] >= 0.1) &
+ (SNRs[iidd] < 0.25))[0]
+ idx3 = np.where((SNRs[iidd] >= 0.25) &
+ (SNRs[iidd] < 0.5))[0]
+ idx4 = np.where((SNRs[iidd] >= 0.5) &
+ (SNRs[iidd] < 0.75))[0]
+ idx5 = np.where((SNRs[iidd] >= 0.75) &
+ (SNRs[iidd] < 0.9))[0]
+ idx6 = np.where(SNRs[iidd] >= 0.9)[0]
+
cv2.drawContours(bkgr_contours, sum([contours[jj] for jj in idx1], []), -1, (255, 0, 0), 1)
cv2.drawContours(bkgr_contours, sum([contours[jj] for jj in idx2], []), -1, (0, 255, 0), 1)
cv2.drawContours(bkgr_contours, sum([contours[jj] for jj in idx3], []), -1, (0, 0, 255), 1)
@@ -105,6 +143,46 @@ def draw_contours():
img.setImage(bkgr_contours, autoLevels=False)
# pg.setConfigOptions(imageAxisOrder='row-major')
+
+
+def draw_contours_update(cf, im):
+ global thrshcomp_line, estimates
+ curFrame = cf.copy()
+
+ if len(estimates.idx_components) > 0:
+ contours = [cv2.findContours(cv2.threshold(img, np.int(thrshcomp_line.value()), 255, 0)[1], cv2.RETR_TREE,
+ cv2.CHAIN_APPROX_SIMPLE)[0] for img in estimates.img_components[estimates.idx_components]]
+ SNRs = np.array(estimates.r_values)
+ iidd = np.array(estimates.idx_components)
+
+ idx1 = np.where(SNRs[iidd] < 0.1)[0]
+ idx2 = np.where((SNRs[iidd] >= 0.1) &
+ (SNRs[iidd] < 0.25))[0]
+ idx3 = np.where((SNRs[iidd] >= 0.25) &
+ (SNRs[iidd] < 0.5))[0]
+ idx4 = np.where((SNRs[iidd] >= 0.5) &
+ (SNRs[iidd] < 0.75))[0]
+ idx5 = np.where((SNRs[iidd] >= 0.75) &
+ (SNRs[iidd] < 0.9))[0]
+ idx6 = np.where(SNRs[iidd] >= 0.9)[0]
+
+ if min_dist_comp in idx1:
+ cv2.drawContours(curFrame, contours[min_dist_comp], -1, (255, 0, 0), 1)
+ if min_dist_comp in idx2:
+ cv2.drawContours(curFrame, contours[min_dist_comp], -1, (0, 255, 0), 1)
+ if min_dist_comp in idx3:
+ cv2.drawContours(curFrame, contours[min_dist_comp], -1, (0, 0, 255), 1)
+ if min_dist_comp in idx4:
+ cv2.drawContours(curFrame, contours[min_dist_comp], -1, (255, 255, 0), 1)
+ if min_dist_comp in idx5:
+ cv2.drawContours(curFrame, contours[min_dist_comp], -1, (255, 0, 255), 1)
+ if min_dist_comp in idx6:
+ cv2.drawContours(curFrame, contours[min_dist_comp], -1, (0, 255, 255), 1)
+
+ im.setImage(curFrame, autoLevels=False)
+
+
+
#%%
@@ -161,6 +239,10 @@ def draw_contours():
layout.addWidget(p2, 1, 1) # plot goes on right side, spanning 2 rows
+#enable only horizontal zoom for the traces component
+p2.setMouseEnabled(x=True, y=False)
+
+
draw_contours()
hist.setLevels(estimates.background_image.min(), estimates.background_image.max())
@@ -181,11 +263,15 @@ def draw_contours():
# zoom to fit imageo
p1.autoRange()
-thrshcomp_line.sigDragged.connect(draw_contours)
+
+mode = "reset"
+p2.setTitle("mode: %s" % (mode))
+
+thrshcomp_line.sigDragged.connect(lambda: draw_contours_overall(mode))
+
def imageHoverEvent(event):
- """Show the position, pixel, and value under the mouse cursor.
- """
+ #Show the position, pixel, and value under the mouse cursor.
global x,y,i,j,val
pos = event.pos()
i, j = pos.y(), pos.x()
@@ -195,32 +281,61 @@ def imageHoverEvent(event):
ppos = img.mapToParent(pos)
x, y = ppos.x(), ppos.y()
-
# Monkey-patch the image to use our custom hover function.
# This is generally discouraged (you should subclass ImageItem instead),
# but it works for a very simple use like this.
img.hoverEvent = imageHoverEvent
+
+
def mouseClickEvent(event):
- global x, y, i, j, val, img2
- distances = np.sum(((x,y)-estimates.cms[estimates.idx_components])**2, axis=1)**0.5
- min_dist_comp = np.argmin(distances)
- estimates.components_to_plot = estimates.idx_components[min_dist_comp]
- p2.plot(estimates.C[estimates.components_to_plot] + estimates.YrA[estimates.components_to_plot], clear=True)
- img2.setImage(estimates.A[:,estimates.components_to_plot].toarray().reshape(estimates.dims, order='F'), autoLevels=True)
- p3.setTitle("pos: (%0.1f, %0.1f) component: %d value: %g dist:%f" % (x, y, estimates.components_to_plot,
- val, distances[min_dist_comp]))
- p1.setTitle("pos: (%0.1f, %0.1f) component: %d value: %g dist:%f" % (x, y, estimates.components_to_plot,
- val, distances[min_dist_comp]))
+ global mode
+ global x,y,i,j,val
+
+ pos = img.mapFromScene(event.pos())
+ x = int(pos.x())
+ y = int(pos.y())
+
+ if x < 0 or x > mov.shape[1] or y < 0 or y > mov.shape[2]:
+ # if the user click outside of the movie, do nothing and jump out of the function
+ return
+
+ i, j = pos.y(), pos.x()
+ i = int(np.clip(i, 0, estimates.background_image.shape[0] - 1))
+ j = int(np.clip(j, 0, estimates.background_image.shape[1] - 1))
+ val = estimates.background_image[i, j, 0]
+
+ if mode is "neurons":
+ show_neurons_clicked()
+
p1.mousePressEvent = mouseClickEvent
+
+#A general rule in Qt is that if you override one mouse event handler, you must override all of them.
+def release(event):
+ pass
+
+p1.mouseReleaseEvent = release
+
+def move(event):
+ pass
+
+p1.mouseMoveEvent = move
+
+
+
+
## 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': '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'}
+ ]
## Create tree of Parameter objects
pars = Parameter.create(name='params', type='group', children=params)
@@ -235,19 +350,146 @@ def mouseClickEvent(event):
{'name': 'REMOVE SINGLE', 'type': 'action'},
{'name': 'SAVE OBJECT', 'type': 'action'}
]
-
+
pars_action = Parameter.create(name='params_action', type='group', children=params_action)
t_action.setParameters(pars_action, showTop=False)
t_action.setWindowTitle('Parameter Action')
+
+def reset_button():
+ global mode
+ mode = "reset"
+ p2.setTitle("mode: %s" % (mode))
+ #clear the upper right image
+ zeros = np.asarray([ [0] * 80 for _ in range(60)])
+ img2.setImage(make_color_img(zeros), autoLevels=False)
+ draw_contours()
+
+
+pars.param('RESET').sigActivated.connect(reset_button)
+
+
+def show_background_button():
+ global bg_vline, min_background, max_background, background_num
+ global mode, background_first_frame_scaled
+ #clear thhe upper right image
+ zeros = np.asarray([ [0] * 80 for _ in range(60)])
+ img2.setImage(make_color_img(zeros), autoLevels=False)
+
+ background_num = (background_num + 1) % estimates.f.shape[0]
+ mode = "background"
+ p2.setTitle("mode: %s %d" % (mode,background_num))
+
+ # display the first frame of the background
+ background_first_frame = estimates.b[:,background_num].reshape(estimates.dims, order='F')
+ min_background_first_frame = np.min(background_first_frame)
+ max_background_first_frame = np.max(background_first_frame)
+ background_first_frame_scaled = make_color_img(background_first_frame, min_max=(min_background_first_frame, max_background_first_frame))
+ img.setImage(background_first_frame_scaled,autoLevels=False)
+
+ # draw the trace and the infinite line
+ trace_background = estimates.f[background_num]
+ p2.plot(trace_background, clear=True)
+ bg_vline = pg.InfiniteLine(angle = 90, movable = True)
+ p2.addItem(bg_vline, ignoreBounds=True)
+ bg_vline.setValue(0)
+ bg_vline.sigPositionChanged.connect(show_background_update)
+
+
+def show_background_update():
+ global bg_index, min_background, max_background, background_scaled
+ bg_index = int(bg_vline.value())
+ if bg_index > -1 and bg_index < estimates.f.shape[-1]:
+ # upper left component scrolls through the frames of the background
+ background = estimates.b[:,background_num].dot(estimates.f[background_num,bg_index]).reshape(estimates.dims, order='F')
+ background_scaled = make_color_img(background, min_max=(min_background[background_num], max_background[background_num]))
+ img.setImage(background_scaled,autoLevels=False)
+
+pars.param('SHOW BACKGROUND').sigActivated.connect(show_background_button)
+
+
+
+
+
+def show_neurons_button():
+ global mode, neuron_selected
+ mode = "neurons"
+ neuron_selected = False
+ p2.setTitle("mode: %s" % (mode))
+ #clear the upper right image
+ zeros = np.asarray([ [0] * 80 for _ in range(60)])
+ img2.setImage(make_color_img(zeros), autoLevels=False)
+
+
+
+def show_neurons_clicked():
+ global nr_vline, nr_index
+ global x,y,i,j,val,min_dist_comp,contour_single, neuron_selected, comp2_scaled
+ neuron_selected = True
+ distances = np.sum(((x,y)-estimates.cms[estimates.idx_components])**2, axis=1)**0.5
+ min_dist_comp = np.argmin(distances)
+ contour_all =[cv2.threshold(img, np.int(thrshcomp_line.value()), 255, 0)[1] for img in estimates.img_components[estimates.idx_components]]
+ contour_single = contour_all[min_dist_comp]
+
+ # draw the traces (lower left component)
+ estimates.components_to_plot = estimates.idx_components[min_dist_comp]
+ p2.plot(estimates.C[estimates.components_to_plot] + estimates.YrA[estimates.components_to_plot], clear=True)
+
+ # plot img (upper left component)
+ img.setImage(estimates.background_image, autoLevels=False)
+ draw_contours_update(estimates.background_image, img)
+ # plot img2 (upper right component)
+ comp2 = np.multiply(estimates.Cn, contour_single > 0)
+ comp2_scaled = make_color_img(comp2, min_max=(np.min(comp2), np.max(comp2)))
+ img2.setImage(comp2_scaled,autoLevels=False)
+ draw_contours_update(comp2_scaled, img2)
+ # set title for the upper two components
+ p3.setTitle("pos: (%0.1f, %0.1f) component: %d value: %g dist:%f" % (x, y, estimates.components_to_plot,
+ val, distances[min_dist_comp]))
+ p1.setTitle("pos: (%0.1f, %0.1f) component: %d value: %g dist:%f" % (x, y, estimates.components_to_plot,
+ val, distances[min_dist_comp]))
+ # draw the infinite line
+ nr_vline = pg.InfiniteLine(angle=90, movable=True)
+ p2.addItem(nr_vline, ignoreBounds=True)
+ nr_vline.setValue(0)
+ nr_vline.sigPositionChanged.connect(show_neurons_update)
+ nr_index = 0
+
+
+def show_neurons_update():
+ global nr_index, frame_denoise_scaled, estimates, raw_mov_scaled
+ global min_mov, max_mov, min_mov_denoise, max_mov_denoise
+ if neuron_selected is False:
+ return
+ nr_index = int(nr_vline.value())
+ if nr_index > 0 and nr_index < mov[:,0,0].shape[0]:
+ # upper left compoenent scrolls through the raw movie
+ raw_mov = mov[nr_index,:,:]
+ raw_mov_scaled = make_color_img(raw_mov, min_max=(min_mov,max_mov))
+ img.setImage(raw_mov_scaled, autoLevels=False)
+ draw_contours_update(raw_mov_scaled, img)
+ # upper right component scrolls through the denoised movie
+ frame_denoise = estimates.A[:,estimates.idx_components].dot(estimates.C[estimates.idx_components,nr_index]).reshape(estimates.dims, order='F')
+ frame_denoise_scaled = make_color_img(frame_denoise, min_max=(min_mov_denoise,max_mov_denoise))
+ img2.setImage(frame_denoise_scaled,autoLevels=False)
+ draw_contours_update(frame_denoise_scaled, img2)
+
+
+pars.param('SHOW NEURONS').sigActivated.connect(show_neurons_button)
+
+
+
+
+
+
+
def add_group():
estimates.accepted_list = np.union1d(estimates.accepted_list,estimates.idx_components)
estimates.rejected_list = np.setdiff1d(estimates.rejected_list,estimates.idx_components)
change(None, None)
-
pars_action.param('ADD GROUP').sigActivated.connect(add_group)
def remove_group():
@@ -290,19 +532,24 @@ def action_pars_activated(param, changes):
## If anything changes in the tree, print a message
def change(param, changes):
- global estimates, cnm_obj, pars, pars_action
- params_obj = cnm_obj.params
+ global estimates, pars, pars_action
set_par = pars.getValues()
if pars_action.param('Filter components').value():
for keyy in set_par.keys():
params_obj.quality.update({keyy: set_par[keyy][0]})
else:
- params_obj.quality.update({'cnn_lowest': 0,'min_cnn_thr':0,'rval_thr':-1,'rval_lowest':-1,'min_SNR':0,'SNR_lowest':0})
-
- estimates.filter_components(mov, params_obj, dview=None,
- select_mode=pars_action.param('View components').value())
-
- draw_contours()
+ params_obj.quality.update({'cnn_lowest': .1,
+ 'min_cnn_thr': 0.99,
+ 'rval_thr': 0.85,
+ 'rval_lowest': -1,
+ 'min_SNR': 2,
+ 'SNR_lowest': 0})
+ estimates.filter_components(mov, params_obj, dview=None,
+ select_mode=pars_action.param('View components').value())
+ if mode is "background":
+ return
+ else:
+ draw_contours()
pars.sigTreeStateChanged.connect(change)
@@ -312,9 +559,6 @@ def change(param, changes):
## END PARAMS
-
-draw_contours()
-
## Display the widget as a new window
w.show()
diff --git a/caiman/mmapping.py b/caiman/mmapping.py
index a4dade990..56a4d41c4 100644
--- a/caiman/mmapping.py
+++ b/caiman/mmapping.py
@@ -1,6 +1,5 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
-
"""
Created on Thu Oct 20 11:33:35 2016
@@ -26,15 +25,17 @@
import caiman as cm
from caiman.paths import memmap_frames_filename
-def prepare_shape(mytuple:Tuple) -> Tuple:
+
+def prepare_shape(mytuple: Tuple) -> Tuple:
""" This promotes the elements inside a shape into np.uint64. It is intended to prevent overflows
with some numpy operations that are sensitive to it, e.g. np.memmap """
if not isinstance(mytuple, tuple):
raise Exception("Internal error: prepare_shape() passed a non-tuple")
return tuple(map(lambda x: np.uint64(x), mytuple))
+
#%%
-def load_memmap(filename:str, mode:str='r') -> Tuple[Any,Tuple,int]:
+def load_memmap(filename: str, mode: str = 'r') -> Tuple[Any, Tuple, int]:
""" Load a memory mapped file created by the function save_memmap
Args:
@@ -63,11 +64,9 @@ def load_memmap(filename:str, mode:str='r') -> Tuple[Any,Tuple,int]:
# 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)
+ 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:
@@ -76,9 +75,19 @@ def load_memmap(filename:str, mode:str='r') -> Tuple[Any,Tuple,int]:
logging.error("Unknown extension for file " + str(filename))
raise Exception('Unknown file extension (should be .mmap)')
+
#%%
-def save_memmap_each(fnames:List[str], dview=None, base_name:str=None, resize_fact=(1, 1, 1), remove_init:int=0,
- idx_xy=None, xy_shifts=None, add_to_movie:float=0, border_to_0:int=0, order:str='C', slices=None) -> List[str]:
+def save_memmap_each(fnames: List[str],
+ dview=None,
+ base_name: str = None,
+ resize_fact=(1, 1, 1),
+ remove_init: int = 0,
+ idx_xy=None,
+ xy_shifts=None,
+ add_to_movie: float = 0,
+ border_to_0: int = 0,
+ order: str = 'C',
+ slices=None) -> List[str]:
"""
Create several memory mapped files using parallel processing
@@ -129,11 +138,16 @@ def save_memmap_each(fnames:List[str], dview=None, base_name:str=None, resize_fa
for idx, f in enumerate(fnames):
if base_name is not None:
- pars.append([f, base_name + '{:04d}'.format(idx), resize_fact[idx], remove_init,
- idx_xy, order, xy_shifts[idx], add_to_movie, border_to_0, slices])
+ pars.append([
+ f, base_name + '{:04d}'.format(idx), resize_fact[idx], remove_init, idx_xy, order, xy_shifts[idx],
+ add_to_movie, border_to_0, slices
+ ])
else:
- pars.append([f, os.path.splitext(f)[0], resize_fact[idx], remove_init, idx_xy, order,
- xy_shifts[idx], add_to_movie, border_to_0, slices])
+ pars.append([
+ f,
+ os.path.splitext(f)[0], resize_fact[idx], remove_init, idx_xy, order, xy_shifts[idx], add_to_movie,
+ border_to_0, slices
+ ])
# Perform the job using whatever computing framework we're set to use
if dview is not None:
@@ -148,7 +162,8 @@ def save_memmap_each(fnames:List[str], dview=None, base_name:str=None, resize_fa
#%%
-def save_memmap_join(mmap_fnames:List[str], base_name:str=None, n_chunks:int=20, dview=None, add_to_mov=0) -> str:
+def save_memmap_join(mmap_fnames: List[str], base_name: str = None, n_chunks: int = 20, dview=None,
+ add_to_mov=0) -> str:
"""
Makes a large file memmap from a number of smaller files
@@ -169,7 +184,7 @@ def save_memmap_join(mmap_fnames:List[str], base_name:str=None, n_chunks:int=20,
order = 'C'
for f in mmap_fnames:
Yr, dims, T = load_memmap(f)
- logging.debug((f, T)) # TODO: Add a text header so this isn't just numeric output, but what to say?
+ logging.debug((f, T)) # TODO: Add a text header so this isn't just numeric output, but what to say?
tot_frames += T
del Yr
@@ -178,15 +193,13 @@ def save_memmap_join(mmap_fnames:List[str], base_name:str=None, n_chunks:int=20,
if base_name is None:
base_name = mmap_fnames[0]
- base_name = base_name[:base_name.find(
- '_d1_')] + '-#-' + str(len(mmap_fnames))
+ base_name = base_name[:base_name.find('_d1_')] + '-#-' + str(len(mmap_fnames))
fname_tot = memmap_frames_filename(base_name, dims, tot_frames, order)
fname_tot = os.path.join(os.path.split(mmap_fnames[0])[0], fname_tot)
logging.info("Memmap file for fname_tot: " + str(fname_tot))
- big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32,
- shape=prepare_shape((d, tot_frames)), order='C')
+ big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32, shape=prepare_shape((d, tot_frames)), order='C')
step = np.int(old_div(d, n_chunks))
pars = []
@@ -194,7 +207,8 @@ def save_memmap_join(mmap_fnames:List[str], base_name:str=None, n_chunks:int=20,
pars.append([fname_tot, d, tot_frames, mmap_fnames, ref, ref + step, add_to_mov])
if len(pars[-1]) != 7:
- raise Exception('You cannot change the number of element in list without changing the statement below (pars[]..)')
+ raise Exception(
+ 'You cannot change the number of element in list without changing the statement below (pars[]..)')
else:
# last batch should include the leftover pixels
pars[-1][-2] = d
@@ -215,6 +229,7 @@ def save_memmap_join(mmap_fnames:List[str], base_name:str=None, n_chunks:int=20,
sys.stdout.flush()
return fname_tot
+
def my_map(dv, func, args) -> List:
v = dv
rc = v.client
@@ -224,7 +239,7 @@ def my_map(dv, func, args) -> List:
amr = v.map(func, args)
pending = set(amr.msg_ids)
- results_all:Dict = dict()
+ results_all: Dict = dict()
counter = 0
while pending:
try:
@@ -236,18 +251,18 @@ def my_map(dv, func, args) -> List:
finished = pending.difference(rc.outstanding)
# update pending to exclude those that just finished
pending = pending.difference(finished)
- if counter%10 == 0:
+ if counter % 10 == 0:
logging.debug(amr.progress)
for msg_id in finished:
# we know these are done, so don't worry about blocking
ar = rc.get_result(msg_id)
- logging.debug("job id %s finished on engine %i" % (msg_id, ar.engine_id))
+ logging.debug("job id %s finished on engine %i" % (msg_id, ar.engine_id)) # TODO: Abstract out the ugly bits
logging.debug("with stdout:")
logging.debug(' ' + ar.stdout.replace('\n', '\n ').rstrip())
logging.debug("and errors:")
logging.debug(' ' + ar.stderr.replace('\n', '\n ').rstrip())
- # note that each job in a map always returns a list of length chunksize
- # even if chunksize == 1
+ # note that each job in a map always returns a list of length chunksize
+ # even if chunksize == 1
results_all.update(ar.get_dict())
counter += 1
@@ -255,25 +270,26 @@ def my_map(dv, func, args) -> List:
del results_all
return result_ordered
+
def save_portion(pars) -> int:
# todo: todocument
use_mmap_save = False
big_mov, d, tot_frames, fnames, idx_start, idx_end, add_to_mov = pars
Ttot = 0
- Yr_tot = np.zeros((idx_end - idx_start, tot_frames), dtype = np.float32)
+ Yr_tot = np.zeros((idx_end - idx_start, tot_frames), dtype=np.float32)
logging.debug("Shape of Yr_tot is " + str(Yr_tot.shape))
for f in fnames:
logging.debug("Saving portion to " + str(f))
Yr, _, T = load_memmap(f)
- Yr_tot[:, Ttot:Ttot + T] = np.ascontiguousarray(Yr[idx_start:idx_end] , dtype = np.float32) + np.float32(add_to_mov)
+ Yr_tot[:, Ttot:Ttot +
+ T] = np.ascontiguousarray(Yr[idx_start:idx_end], dtype=np.float32) + np.float32(add_to_mov)
Ttot = Ttot + T
del Yr
logging.debug("Index start and end are " + str(idx_start) + " and " + str(idx_end))
if use_mmap_save:
- big_mov = np.memmap(big_mov, mode='r+', dtype=np.float32,
- shape=prepare_shape((d, tot_frames)), order='C')
+ big_mov = np.memmap(big_mov, mode='r+', dtype=np.float32, shape=prepare_shape((d, tot_frames)), order='C')
big_mov[idx_start:idx_end, :] = Yr_tot
del big_mov
else:
@@ -282,10 +298,12 @@ def save_portion(pars) -> int:
f.write(Yr_tot)
computed_position = np.uint64(idx_end * Yr_tot.dtype.itemsize * tot_frames)
if f.tell() != computed_position:
- logging.critical(f"Error in mmap portion write: at position {f.tell()}")
- logging.critical(f"But should be at position {idx_end} * {Yr_tot.dtype.itemsize} * {tot_frames} = {computed_position}")
- f.close()
- raise Exception('Internal error in mmapping: Actual position does not match computed position')
+ logging.critical(f"Error in mmap portion write: at position {f.tell()}")
+ logging.critical(
+ f"But should be at position {idx_end} * {Yr_tot.dtype.itemsize} * {tot_frames} = {computed_position}"
+ )
+ f.close()
+ raise Exception('Internal error in mmapping: Actual position does not match computed position')
del Yr_tot
logging.debug('done')
@@ -293,24 +311,39 @@ def save_portion(pars) -> int:
#%%
-def save_place_holder(pars:List) -> str:
+def save_place_holder(pars: List) -> str:
""" To use map reduce
"""
# todo: todocument
- (f, base_name, resize_fact, remove_init, idx_xy, order,
- xy_shifts, add_to_movie, border_to_0, slices) = pars
+ (f, base_name, resize_fact, remove_init, idx_xy, order, xy_shifts, add_to_movie, border_to_0, slices) = pars
- return save_memmap([f], base_name=base_name, resize_fact=resize_fact, remove_init=remove_init,
- idx_xy=idx_xy, order=order, xy_shifts=xy_shifts,
- add_to_movie=add_to_movie, border_to_0=border_to_0, slices=slices)
+ return save_memmap([f],
+ base_name=base_name,
+ resize_fact=resize_fact,
+ remove_init=remove_init,
+ idx_xy=idx_xy,
+ order=order,
+ xy_shifts=xy_shifts,
+ add_to_movie=add_to_movie,
+ border_to_0=border_to_0,
+ slices=slices)
#%%
-def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1, 1), remove_init:int=0, idx_xy:Tuple=None,
- order:str='F', xy_shifts:Optional[List]=None, is_3D:bool=False, add_to_movie:float=0, border_to_0=0, dview = None,
- n_chunks:int=100, slices=None) -> str:
-
+def save_memmap(filenames: List[str],
+ base_name: str = 'Yr',
+ resize_fact: Tuple = (1, 1, 1),
+ remove_init: int = 0,
+ idx_xy: Tuple = None,
+ order: str = 'F',
+ xy_shifts: Optional[List] = None,
+ is_3D: bool = False,
+ add_to_movie: float = 0,
+ border_to_0=0,
+ dview=None,
+ n_chunks: int = 100,
+ slices=None) -> str:
""" Efficiently write data from a list of tif files into a memory mappable file
Args:
@@ -369,7 +402,7 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1
recompute_each_memmap = False
for file__ in filenames:
if ('order_' + order not in file__) or ('.mmap' not in file__):
- recompute_each_memmap = True
+ recompute_each_memmap = True
if recompute_each_memmap or (remove_init>0) or (idx_xy is not None)\
@@ -379,16 +412,16 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1
logging.debug('Distributing memory map over many files')
# Here we make a bunch of memmap files in the right order. Same parameters
fname_parts = cm.save_memmap_each(filenames,
- base_name = base_name,
- order = order,
- border_to_0 = border_to_0,
- dview = dview,
- resize_fact = resize_fact,
- remove_init = remove_init,
- idx_xy = idx_xy,
- xy_shifts = xy_shifts,
- slices = slices,
- add_to_movie = add_to_movie)
+ base_name=base_name,
+ order=order,
+ border_to_0=border_to_0,
+ dview=dview,
+ resize_fact=resize_fact,
+ remove_init=remove_init,
+ idx_xy=idx_xy,
+ xy_shifts=xy_shifts,
+ slices=slices,
+ add_to_movie=add_to_movie)
else:
fname_parts = filenames
@@ -396,26 +429,25 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1
if order == 'F':
raise Exception('You cannot merge files in F order, they must be in C order for CaImAn')
-
fname_new = cm.save_memmap_join(fname_parts, base_name=base_name, dview=dview, n_chunks=n_chunks)
else:
- # TODO: can be done online
+ # TODO: can be done online
Ttot = 0
for idx, f in enumerate(filenames):
- if isinstance(f, str): # Might not always be filenames.
+ if isinstance(f, str): # Might not always be filenames.
logging.debug(f)
if is_3D:
- Yr = f if not(isinstance(f, basestring)) else tifffile.imread(f)
+ Yr = f if not (isinstance(f, basestring)) else tifffile.imread(f)
if slices is not None:
Yr = Yr[slices]
else:
- if idx_xy is None: #todo remove if not used, superceded by the slices parameter
+ if idx_xy is None: #todo remove if not used, superceded by the slices parameter
Yr = Yr[remove_init:]
- elif len(idx_xy) == 2: #todo remove if not used, superceded by the slices parameter
+ elif len(idx_xy) == 2: #todo remove if not used, superceded by the slices parameter
Yr = Yr[remove_init:, idx_xy[0], idx_xy[1]]
- else: #todo remove if not used, superceded by the slices parameter
+ else: #todo remove if not used, superceded by the slices parameter
Yr = Yr[remove_init:, idx_xy[0], idx_xy[1], idx_xy[2]]
else:
@@ -439,9 +471,11 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1
Yr = np.array(Yr)[remove_init:, idx_xy[0], idx_xy[1], idx_xy[2]]
if border_to_0 > 0:
- if slices is not None:
+ if slices is not None:
if type(slices) is list:
- raise Exception('You cannot slice in x and y and then use add_to_movie: if you only want to slice in time do not pass in a list but just a slice object')
+ raise Exception(
+ 'You cannot slice in x and y and then use add_to_movie: if you only want to slice in time do not pass in a list but just a slice object'
+ )
min_mov = Yr.calc_min()
Yr[:, :border_to_0, :] = min_mov
@@ -461,21 +495,28 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1
Yr = np.ascontiguousarray(Yr, dtype=np.float32) + np.float32(0.0001) + np.float32(add_to_movie)
if idx == 0:
- fname_tot = base_name + '_d1_' + str(dims[0]) + '_d2_' + str(dims[1]) + '_d3_' + str(
- 1 if len(dims) == 2 else dims[2]) + '_order_' + str(order) # TODO: Rewrite more legibly
+ fname_tot = base_name + '_d1_' + str(
+ dims[0]) + '_d2_' + str(dims[1]) + '_d3_' + str(1 if len(dims) == 2 else dims[2]) + '_order_' + str(
+ order) # TODO: Rewrite more legibly
if isinstance(f, str):
fname_tot = os.path.join(os.path.split(f)[0], fname_tot)
if len(filenames) > 1:
- big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32,
- shape=prepare_shape((np.prod(dims), T)), order=order)
+ big_mov = np.memmap(fname_tot,
+ mode='w+',
+ dtype=np.float32,
+ shape=prepare_shape((np.prod(dims), T)),
+ order=order)
big_mov[:, Ttot:Ttot + T] = Yr
del big_mov
else:
logging.debug('SAVING WITH numpy.tofile()')
Yr.tofile(fname_tot)
else:
- big_mov = np.memmap(fname_tot, dtype=np.float32, mode='r+',
- shape=prepare_shape((np.prod(dims), Ttot + T)), order=order)
+ big_mov = np.memmap(fname_tot,
+ dtype=np.float32,
+ mode='r+',
+ shape=prepare_shape((np.prod(dims), Ttot + T)),
+ order=order)
big_mov[:, Ttot:Ttot + T] = Yr
del big_mov
@@ -493,10 +534,12 @@ def save_memmap(filenames:List[str], base_name:str='Yr', resize_fact:Tuple=(1, 1
return fname_new
+
#%%
-def parallel_dot_product(A:np.ndarray, b, block_size:int=5000, dview=None, transpose=False, num_blocks_per_run=20) -> np.ndarray:
+def parallel_dot_product(A: np.ndarray, b, block_size: int = 5000, dview=None, transpose=False,
+ num_blocks_per_run=20) -> np.ndarray:
# todo: todocument
""" Chunk matrix product between matrix and column vectors
@@ -550,11 +593,9 @@ def parallel_dot_product(A:np.ndarray, b, block_size:int=5000, dview=None, trans
for itera in range(0, len(pars), num_blocks_per_run):
if 'multiprocessing' in str(type(dview)):
- results = dview.map_async(
- dot_place_holder, pars[itera:itera + num_blocks_per_run]).get(4294967)
+ results = dview.map_async(dot_place_holder, pars[itera:itera + num_blocks_per_run]).get(4294967)
else:
- results = dview.map_sync(
- dot_place_holder, pars[itera:itera + num_blocks_per_run])
+ results = dview.map_sync(dot_place_holder, pars[itera:itera + num_blocks_per_run])
logging.debug('Processed:' + str([itera, itera + len(results)]))
@@ -566,7 +607,6 @@ def parallel_dot_product(A:np.ndarray, b, block_size:int=5000, dview=None, trans
else:
logging.debug('Filling')
-
for res in results:
output[res[0]] = res[1]
@@ -577,7 +617,7 @@ def parallel_dot_product(A:np.ndarray, b, block_size:int=5000, dview=None, trans
#%%
-def dot_place_holder(par:List) -> Tuple:
+def dot_place_holder(par: List) -> Tuple:
# todo: todocument
A_name, idx_to_pass, b_, transpose = par
@@ -587,9 +627,9 @@ def dot_place_holder(par:List) -> Tuple:
logging.debug((idx_to_pass[-1]))
if 'sparse' in str(type(b_)):
if transpose:
- outp = (b_.T.tocsc()[:, idx_to_pass].dot(
- A_[idx_to_pass])).T.astype(np.float32)
-
+ # outp = (b_.tocsr()[idx_to_pass].T.dot(
+ # A_[idx_to_pass])).T.astype(np.float32)
+ outp = (b_.T.tocsc()[:, idx_to_pass].dot(A_[idx_to_pass])).T.astype(np.float32)
else:
outp = (b_.T.dot(A_[idx_to_pass].T)).T.astype(np.float32)
else:
@@ -603,23 +643,25 @@ def dot_place_holder(par:List) -> Tuple:
#%%
-def save_tif_to_mmap_online(movie_iterable, save_base_name='YrOL_', order='C',
- add_to_movie=0, border_to_0=0) -> str:
+def save_tif_to_mmap_online(movie_iterable, save_base_name='YrOL_', order='C', add_to_movie=0, border_to_0=0) -> str:
# todo: todocument
- if isinstance(movie_iterable, basestring): # Allow specifying a filename rather than its data rep
- with tifffile.TiffFile(movie_iterable) as tf: # And load it if that happens
+ if isinstance(movie_iterable, basestring): # Allow specifying a filename rather than its data rep
+ with tifffile.TiffFile(movie_iterable) as tf: # And load it if that happens
movie_iterable = cm.movie(tf)
count = 0
new_mov = []
- dims = (len(movie_iterable),) + movie_iterable[0].shape # TODO: Don't pack length into dims
+ dims = (len(movie_iterable),) + movie_iterable[0].shape # TODO: Don't pack length into dims
fname_tot = (memmap_frames_filename(save_base_name, dims[1:], dims[0], order))
- big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32,
- shape=prepare_shape((np.prod(dims[1:]), dims[0])), order=order)
+ big_mov = np.memmap(fname_tot,
+ mode='w+',
+ dtype=np.float32,
+ shape=prepare_shape((np.prod(dims[1:]), dims[0])),
+ order=order)
for page in movie_iterable:
if count % 100 == 0:
@@ -631,8 +673,7 @@ def save_tif_to_mmap_online(movie_iterable, save_base_name='YrOL_', order='C',
img = np.array(page, dtype=np.float32)
new_img = img
if save_base_name is not None:
- big_mov[:, count] = np.reshape(
- new_img, np.prod(dims[1:]), order='F')
+ big_mov[:, count] = np.reshape(new_img, np.prod(dims[1:]), order='F')
else:
new_mov.append(new_img)
@@ -642,8 +683,7 @@ def save_tif_to_mmap_online(movie_iterable, save_base_name='YrOL_', order='C',
img[:, -border_to_0:] = 0
img[-border_to_0:, :] = 0
- big_mov[:, count] = np.reshape(
- img + add_to_movie, np.prod(dims[1:]), order='F')
+ big_mov[:, count] = np.reshape(img + add_to_movie, np.prod(dims[1:]), order='F')
count += 1
big_mov.flush()
diff --git a/caiman/motion_correction.py b/caiman/motion_correction.py
index 75e533ab1..15f12fd49 100644
--- a/caiman/motion_correction.py
+++ b/caiman/motion_correction.py
@@ -215,7 +215,7 @@ def motion_correct(self, template=None, save_movie=False):
output can be saved as a memory mapped file.
Args:
- template: nd.array, default: None
+ template: ndarray, default: None
template provided by user for motion correction
save_movie: bool, default: False
@@ -253,7 +253,7 @@ def motion_correct(self, template=None, save_movie=False):
self.mmap_file = self.fname_tot_els if self.pw_rigid else self.fname_tot_rig
return self
- def motion_correct_rigid(self, template=None, save_movie=False):
+ def motion_correct_rigid(self, template=None, save_movie=False) -> None:
"""
Perform rigid motion correction
@@ -264,9 +264,6 @@ def motion_correct_rigid(self, template=None, save_movie=False):
save_movie_rigid:Bool
save the movies vs just get the template
- Returns:
- self
-
Important Fields:
self.fname_tot_rig: name of the mmap file saved
@@ -308,28 +305,19 @@ def motion_correct_rigid(self, template=None, save_movie=False):
self.fname_tot_rig += [_fname_tot_rig]
self.shifts_rig += _shifts_rig
- return self
-
- def motion_correct_pwrigid(
- self,
- save_movie=True,
- template=None,
- show_template=False):
+ def motion_correct_pwrigid(self, save_movie:bool=True, template:np.ndarray=None, show_template:bool=False) -> None:
"""Perform pw-rigid motion correction
Args:
- template: ndarray 2D (or 3D)
- if known, one can pass a template to register the frames to
-
save_movie:Bool
save the movies vs just get the template
+ template: ndarray 2D (or 3D)
+ if known, one can pass a template to register the frames to
+
show_template: boolean
whether to show the updated template at each iteration
- Returns:
- self
-
Important Fields:
self.fname_tot_els: name of the mmap file saved
self.templates_els: template updated by iterating over the chunks
@@ -347,7 +335,7 @@ def motion_correct_pwrigid(
num_iter = 1
if template is None:
logging.info('Generating template by rigid motion correction')
- self = self.motion_correct_rigid()
+ self.motion_correct_rigid()
self.total_template_els = self.total_template_rig.copy()
else:
self.total_template_els = template
@@ -388,10 +376,9 @@ def motion_correct_pwrigid(
if self.is3D:
self.z_shifts_els += _z_shifts_els
self.coord_shifts_els += _coord_shifts_els
-# return self
- def apply_shifts_movie(self, fname, rigid_shifts=True, save_memmap=False,
- save_base_name='MC', order='F'):
+ def apply_shifts_movie(self, fname, rigid_shifts:bool=True, save_memmap:bool=False,
+ save_base_name:str='MC', order:str='F'):
"""
Applies shifts found by registering one file to a different file. Useful
for cases when shifts computed from a structural channel are applied to a
@@ -443,25 +430,26 @@ def apply_shifts_movie(self, fname, rigid_shifts=True, save_memmap=False,
# np.float32), np.arange(0., dims[1]).astype(np.float32))
x_grid, y_grid = np.meshgrid(np.arange(0., dims[1]).astype(
np.float32), np.arange(0., dims[0]).astype(np.float32))
- m_reg = [cv2.remap(img,
- -cv2.resize(shiftY, dims[::-1]) + x_grid, -cv2.resize(shiftX, dims[::-1]) + y_grid, cv2.INTER_CUBIC)
+ m_reg = [cv2.remap(img, -cv2.resize(shiftY, dims[::-1]) + x_grid,
+ -cv2.resize(shiftX, dims[::-1]) + y_grid,
+ cv2.INTER_CUBIC, borderMode=cv2.BORDER_REPLICATE)
for img, shiftX, shiftY in zip(Y, shifts_x, shifts_y)]
- if save_memmap:
- Y = np.array(Y)
- dims = Y.shape
- fname_tot = memmap_frames_filename(save_base_name, dims[1:], dims[0], order)
- big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32,
- shape=prepare_shape((np.prod(dims[1:]), dims[0])), order=order)
- big_mov[:] = np.reshape(Y.transpose(1, 2, 0), (np.prod(dims[1:]), dims[0]), order=order)
- big_mov.flush()
- del big_mov
- return fname_tot
- else:
- return cm.movie(np.stack(m_reg, axis=0))
+ m_reg = np.stack(m_reg, axis=0)
+ if save_memmap:
+ dims = m_reg.shape
+ fname_tot = memmap_frames_filename(save_base_name, dims[1:], dims[0], order)
+ big_mov = np.memmap(fname_tot, mode='w+', dtype=np.float32,
+ shape=prepare_shape((np.prod(dims[1:]), dims[0])), order=order)
+ big_mov[:] = np.reshape(m_reg.transpose(1, 2, 0), (np.prod(dims[1:]), dims[0]), order=order)
+ big_mov.flush()
+ del big_mov
+ return fname_tot
+ else:
+ return cm.movie(m_reg)
#%%
-def apply_shift_iteration(img, shift, border_nan=False, border_type=cv2.BORDER_REFLECT):
+def apply_shift_iteration(img, shift, border_nan:bool=False, border_type=cv2.BORDER_REFLECT):
# todo todocument
sh_x_n, sh_y_n = shift
@@ -1304,8 +1292,8 @@ def init_cuda_process():
cudadrv.init()
dev = cudadrv.Device(0)
- cudactx = dev.make_context()
- atexit.register(cudactx.pop)
+ cudactx = dev.make_context() # type: ignore
+ atexit.register(cudactx.pop) # type: ignore
def close_cuda_process(n):
@@ -1313,9 +1301,11 @@ def close_cuda_process(n):
Cleanup cuda process
"""
+ global cudactx
+
import skcuda.misc as cudamisc
try:
- cudamisc.done_context(cudactx)
+ cudamisc.done_context(cudactx) # type: ignore
except:
pass
@@ -1586,7 +1576,7 @@ def register_translation(src_image, target_image, upsample_factor=1,
from skcuda.fft import fft as cudafft
from skcuda.fft import ifft as cudaifft
try:
- cudactx
+ cudactx # type: ignore
except NameError:
init_cuda_process()
diff --git a/caiman/paths.py b/caiman/paths.py
index 92b2de63b..d50cb3360 100644
--- a/caiman/paths.py
+++ b/caiman/paths.py
@@ -1,5 +1,4 @@
#!/usr/bin/env python
-
""" Utilities for retrieving paths around the caiman package and its datadirs
"""
@@ -7,22 +6,25 @@
import os
from typing import Tuple
+
#######
# datadir
def caiman_datadir() -> str:
- """
+ """
The datadir is a user-configurable place which holds a user-modifiable copy of
data the Caiman libraries need to function, alongside code demos and other things.
This is meant to be separate from the library install of Caiman, which may be installed
into the global python library path (or into a conda path or somewhere else messy).
"""
- if "CAIMAN_DATA" in os.environ:
- return os.environ["CAIMAN_DATA"]
- else:
- return os.path.join(os.path.expanduser("~"), "caiman_data")
+ if "CAIMAN_DATA" in os.environ:
+ return os.environ["CAIMAN_DATA"]
+ else:
+ return os.path.join(os.path.expanduser("~"), "caiman_data")
+
def caiman_datadir_exists() -> bool:
- return os.path.isdir(caiman_datadir())
+ return os.path.isdir(caiman_datadir())
+
######
# memmap files
@@ -30,14 +32,14 @@ def caiman_datadir_exists() -> bool:
# Right now these are usually stored in the cwd of the script, although the basename could change that
# In the future we may consistently store these somewhere under the caiman_datadir
-def memmap_frames_filename(basename:str, dims:Tuple, frames:int, order:str='F') -> str:
- # Some functions calling this have the first part of *their* dims Tuple be the number of frames.
- # They *must* pass a slice to this so dims is only X, Y, and optionally Z. Frames is passed separately.
- dimfield_0 = dims[0]
- dimfield_1 = dims[1]
- if len(dims) == 3:
- dimfield_2 = dims[2]
- else:
- dimfield_2 = 1
- return f"{basename}_d1_{dimfield_0}_d2_{dimfield_1}_d3_{dimfield_2}_order_{order}_frames_{frames}_.mmap"
+def memmap_frames_filename(basename: str, dims: Tuple, frames: int, order: str = 'F') -> str:
+ # Some functions calling this have the first part of *their* dims Tuple be the number of frames.
+ # They *must* pass a slice to this so dims is only X, Y, and optionally Z. Frames is passed separately.
+ dimfield_0 = dims[0]
+ dimfield_1 = dims[1]
+ if len(dims) == 3:
+ dimfield_2 = dims[2]
+ else:
+ dimfield_2 = 1
+ return f"{basename}_d1_{dimfield_0}_d2_{dimfield_1}_d3_{dimfield_2}_order_{order}_frames_{frames}_.mmap"
diff --git a/caiman/source_extraction/cnmf/cnmf.py b/caiman/source_extraction/cnmf/cnmf.py
index e7b9dd95a..84eb845fd 100644
--- a/caiman/source_extraction/cnmf/cnmf.py
+++ b/caiman/source_extraction/cnmf/cnmf.py
@@ -76,7 +76,7 @@ class CNMF(object):
"""
def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p=2, dview=None,
Ain=None, Cin=None, b_in=None, f_in=None, do_merge=True,
- ssub=2, tsub=2, p_ssub=1, p_tsub=1, method_init='greedy_roi', alpha_snmf=None,
+ ssub=2, tsub=2, p_ssub=1, p_tsub=1, method_init='greedy_roi', alpha_snmf=100,
rf=None, stride=None, memory_fact=1, gnb=1, nb_patch=1, only_init_patch=False,
method_deconvolution='oasis', n_pixels_per_process=4000, block_size_temp=5000, num_blocks_per_run_temp=20,
block_size_spat=5000, num_blocks_per_run_spat=20,
@@ -94,7 +94,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p
max_num_added=3, min_num_trial=2, thresh_CNN_noisy=0.5,
fr=30, decay_time=0.4, min_SNR=2.5, ssub_B=2, init_iter=2,
sniper_mode=False, use_peak_max=False, test_both=False,
- expected_comps=500, params=None):
+ expected_comps=500, max_merge_area=None, params=None):
"""
Constructor of the CNMF method
@@ -112,7 +112,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p
merging threshold, max correlation allowed
dview: Direct View object
- for parallelization pruposes when using ipyparallel
+ for parallelization purposes when using ipyparallel
p: int
order of the autoregressive process used to estimate deconvolution
@@ -193,7 +193,7 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p
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.
- Howeverm benefits can be considerable if done because if many components (>2000) are created
+ 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
@@ -254,6 +254,9 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p
init_iter: int, optional
number of iterations for 1-photon imaging initialization
+
+ max_merge_area: int, optional
+ maximum area (in pixels) of merged components, used to determine whether to merge components during fitting process
"""
self.dview = dview
@@ -289,7 +292,8 @@ def __init__(self, n_processes, k=5, gSig=[4, 4], gSiz=None, merge_thresh=0.8, p
n_refit=n_refit, num_times_comp_updated=num_times_comp_updated, simultaneously=simultaneously,
sniper_mode=sniper_mode, test_both=test_both, thresh_CNN_noisy=thresh_CNN_noisy,
thresh_fitness_delta=thresh_fitness_delta, thresh_fitness_raw=thresh_fitness_raw, thresh_overlap=thresh_overlap,
- update_num_comps=update_num_comps, use_dense=use_dense, use_peak_max=use_peak_max
+ update_num_comps=update_num_comps, use_dense=use_dense, use_peak_max=use_peak_max, alpha_snmf=alpha_snmf,
+ max_merge_area=max_merge_area
)
else:
self.params = params
@@ -504,7 +508,7 @@ def fit(self, images, indices=[slice(None), slice(None)]):
self.params.set('temporal', {'p': 0})
else:
self.params.set('temporal', {'p': self.params.get('preprocess', 'p')})
- logging.info('deconvolution ...')
+ logging.info('deconvolution ...')
self.update_temporal(Yr)
@@ -512,7 +516,7 @@ def fit(self, images, indices=[slice(None), slice(None)]):
logging.info('refinement...')
if self.params.get('merging', 'do_merge'):
logging.info('merging components ...')
- self.merge_comps(Yr, mx=50, fast_merge=True)
+ self.merge_comps(Yr, mx=50, fast_merge=True, max_merge_area=self.params.get('merging', 'max_merge_area'))
logging.info('Updating spatial ...')
@@ -584,11 +588,11 @@ def fit(self, images, indices=[slice(None), slice(None)]):
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]])])
+ #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]])])
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,
@@ -676,6 +680,7 @@ def compute_residuals(self, Yr):
AA = Ab.T.dot(Ab) * nA2_inv_mat
self.estimates.YrA = (YA - (AA.T.dot(Cf)).T)[:, :self.estimates.A.shape[-1]].T
+ self.estimates.R = self.estimates.YrA
return self
@@ -855,6 +860,7 @@ def update_temporal(self, Y, use_init=True, **kwargs):
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):
@@ -890,17 +896,19 @@ def update_spatial(self, Y, use_init=True, **kwargs):
return self
- def merge_comps(self, Y, mx=50, fast_merge=True):
+ def merge_comps(self, Y, mx=50, fast_merge=True, max_merge_area=None):
"""merges components
"""
self.estimates.A, self.estimates.C, self.estimates.nr, self.estimates.merged_ROIs, self.estimates.S, \
- self.estimates.bl, self.estimates.c1, self.estimates.neurons_sn, self.estimates.g, self.empty_merged=\
- merge_components(Y, self.estimates.A, self.estimates.b, self.estimates.C, self.estimates.f, self.estimates.S,
- self.estimates.sn, self.params.get_group('temporal'),
+ self.estimates.bl, self.estimates.c1, self.estimates.neurons_sn, self.estimates.g, self.empty_merged, \
+ self.estimates.YrA =\
+ merge_components(Y, self.estimates.A, self.estimates.b, self.estimates.C, self.estimates.YrA,
+ self.estimates.f, self.estimates.S, self.estimates.sn, self.params.get_group('temporal'),
self.params.get_group('spatial'), dview=self.dview,
bl=self.estimates.bl, c1=self.estimates.c1, sn=self.estimates.neurons_sn,
g=self.estimates.g, thr=self.params.get('merging', 'merge_thr'), mx=mx,
- fast_merge=fast_merge)
+ fast_merge=fast_merge, merge_parallel=self.params.get('merging', 'merge_parallel'),
+ max_merge_area=max_merge_area)
return self
@@ -950,41 +958,101 @@ def load_CNMF(filename, n_processes=1, dview=None):
Args:
filename: str
- hdf5 file name containing the saved object
- dview: multiprocessingor ipyparallel object
+ hdf5 (or nwb) file name containing the saved object
+ dview: multiprocessing or ipyparallel object
useful to set up parllelization in the objects
'''
new_obj = CNMF(n_processes)
- for key, val in load_dict_from_hdf5(filename).items():
- if key == 'params':
- prms = CNMFParams()
- prms.spatial = val['spatial']
- prms.temporal = val['temporal']
- prms.patch = val['patch']
- prms.preprocess = val['preprocess']
- prms.init = val['init']
- prms.merging = val['merging']
- prms.quality = val['quality']
- prms.data = val['data']
- prms.online = val['online']
- prms.motion = val['motion']
- setattr(new_obj, key, prms)
- elif key == 'dview':
- setattr(new_obj, key, dview)
- elif key == 'estimates':
- estims = Estimates()
- for kk, vv in val.items():
- if kk == 'discarded_components':
- if vv is not None:
- discarded_components = Estimates()
- for kk__, vv__ in vv.items():
- setattr(discarded_components, kk__, vv__)
- setattr(estims, kk, discarded_components)
- else:
- setattr(estims, kk, vv)
+ if os.path.splitext(filename)[1].lower() in ('.hdf5', '.h5'):
+ for key, val in load_dict_from_hdf5(filename).items():
+ if key == 'params':
+ prms = CNMFParams()
+ prms.spatial = val['spatial']
+ prms.temporal = val['temporal']
+ prms.patch = val['patch']
+ prms.preprocess = val['preprocess']
+ prms.init = val['init']
+ prms.merging = val['merging']
+ prms.quality = val['quality']
+ prms.data = val['data']
+ prms.online = val['online']
+ prms.motion = val['motion']
+ setattr(new_obj, key, prms)
+ elif key == 'dview':
+ setattr(new_obj, key, dview)
+ elif key == 'estimates':
+ estims = Estimates()
+ for kk, vv in val.items():
+ if kk == 'discarded_components':
+ if vv is not None:
+ discarded_components = Estimates()
+ for kk__, vv__ in vv.items():
+ setattr(discarded_components, kk__, vv__)
+ setattr(estims, kk, discarded_components)
+ else:
+ setattr(estims, kk, vv)
- setattr(new_obj, key, estims)
- else:
- setattr(new_obj, key, val)
+ setattr(new_obj, key, estims)
+ else:
+ setattr(new_obj, key, val)
+ elif os.path.splitext(filename)[1].lower() == '.nwb':
+ from pynwb import NWBHDF5IO
+ with NWBHDF5IO(filename, 'r') as io:
+ nwb = io.read()
+ ophys = nwb.processing['ophys']
+ rrs_group = ophys.data_interfaces['Fluorescence'].roi_response_series
+ rrs = rrs_group['RoiResponseSeries']
+ C = rrs.data[:].T
+ rois = rrs.rois
+ roi_indices = rois.data
+ A = rois.table['image_mask'][roi_indices, ...]
+ dims = A.shape[1:]
+ A = A.reshape((A.shape[0], -1)).T
+ A = scipy.sparse.csc_matrix(A)
+ if 'Background_Fluorescence_Response' in rrs_group:
+ brs = rrs_group['Background_Fluorescence_Response']
+ f = brs.data[:].T
+ brois = brs.rois
+ broi_indices = brois.data
+ b = brois.table['image_mask'][broi_indices, ...]
+ b = b.reshape((b.shape[0], -1)).T
+ 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
+
+ frame_rate = ophys.data_interfaces['ImageSegmentation'].plane_segmentations['PlaneSegmentation']. \
+ imaging_plane.imaging_rate
+
+ if 'r' in rois.table:
+ estims.r_values = rois.table['r'][roi_indices]
+ if 'snr' in rois.table:
+ estims.SNR_comp = rois.table['snr'][roi_indices]
+ if 'cnn' in rois.table:
+ estims.cnn_preds = rois.table['cnn'][roi_indices]
+ if 'keep' in rois.table:
+ keep = rois.table['keep'][roi_indices]
+ estims.idx_components = np.where(keep)[0]
+ estims.nr = len(roi_indices)
+
+ if 'summary_images' in ophys.data_interfaces:
+ if 'Cn' in ophys.data_interfaces['summary_images']:
+ estims.Cn = ophys.data_interfaces['summary_images']['Cn']
+ if hasattr(nwb.acquisition['TwoPhotonSeries'], 'external_file'):
+ setattr(new_obj, 'mmap_file', nwb.acquisition['TwoPhotonSeries'].external_file[0])
+ else:
+ setattr(new_obj, 'mmap_file', filename)
+
+ estims.dims = dims
+ prms = CNMFParams(dims=dims)
+ prms.set('data', {'fr': frame_rate})
+
+ setattr(new_obj, 'params', prms)
+ setattr(new_obj, 'dview', dview)
+ setattr(new_obj, 'estimates', estims)
+
+ else:
+ raise NotImplementedError('unsupported file extension')
return new_obj
diff --git a/caiman/source_extraction/cnmf/estimates.py b/caiman/source_extraction/cnmf/estimates.py
index 86c905bce..f4a0f7d25 100644
--- a/caiman/source_extraction/cnmf/estimates.py
+++ b/caiman/source_extraction/cnmf/estimates.py
@@ -11,15 +11,19 @@
import numpy as np
import scipy.sparse
from typing import List
+import time
import caiman
from .utilities import detrend_df_f
from .spatial import threshold_components
+from .temporal import constrained_foopsi_parallel
from .merging import merge_iteration
from ...components_evaluation import (
evaluate_components_CNN, estimate_components_quality_auto,
select_components_from_metrics)
-from ...base.rois import detect_duplicates_and_subsets, nf_match_neurons_in_binary_masks, nf_masks_to_neurof_dict
+from ...base.rois import (
+ detect_duplicates_and_subsets, nf_match_neurons_in_binary_masks,
+ nf_masks_to_neurof_dict)
from .initialization import downscale
@@ -180,7 +184,7 @@ def plot_contours(self, img=None, idx=None, crd=None, thr_method='max',
if img is None:
img = np.reshape(np.array(self.A.mean(1)), self.dims, order='F')
if self.coordinates is None: # not hasattr(self, 'coordinates'):
- self.coordinates = caiman.utils.visualization.get_contours(self.A, self.dims, thr=thr, thr_method=thr_method)
+ self.coordinates = caiman.utils.visualization.get_contours(self.A, img.shape, thr=thr, thr_method=thr_method)
plt.figure()
if params is not None:
plt.suptitle('min_SNR=%1.2f, rval_thr=%1.2f, use_cnn=%i'
@@ -378,6 +382,55 @@ def nb_view_components(self, Yr=None, img=None, idx=None,
thr=thr, denoised_color=denoised_color, cmap=cmap)
return self
+ def hv_view_components(self, Yr=None, img=None, idx=None,
+ denoised_color=None, cmap='viridis'):
+ """view spatial and temporal components interactively in a notebook
+
+ Args:
+ Yr : np.ndarray
+ movie in format pixels (d) x frames (T)
+
+ img : np.ndarray
+ background image for contour plotting. Default is the mean
+ image of all spatial components (d1 x d2)
+
+ idx : list
+ list of components to be plotted
+
+ denoised_color: string or None
+ color name (e.g. 'red') or hex color code (e.g. '#F0027F')
+
+ cmap: string
+ name of colormap (e.g. 'viridis') used to plot image_neurons
+ """
+ if 'csc_matrix' not in str(type(self.A)):
+ self.A = scipy.sparse.csc_matrix(self.A)
+
+ plt.ion()
+ nr, T = self.C.shape
+ if self.R is None:
+ self.R = self.YrA
+ if self.R.shape != [nr, T]:
+ if self.YrA is None:
+ self.compute_residuals(Yr)
+ else:
+ self.R = self.YrA
+
+ if img is None:
+ img = np.reshape(np.array(self.A.mean(axis=1)), self.dims, order='F')
+
+ if idx is None:
+ hv_plot = caiman.utils.visualization.hv_view_patches(
+ Yr, self.A, self.C, self.b, self.f, self.dims[0], self.dims[1],
+ YrA=self.R, image_neurons=img, denoised_color=denoised_color,
+ cmap=cmap)
+ else:
+ hv_plot = caiman.utils.visualization.hv_view_patches(
+ Yr, self.A.tocsc()[:, idx], self.C[idx], self.b, self.f,
+ self.dims[0], self.dims[1], YrA=self.R[idx], image_neurons=img,
+ denoised_color=denoised_color, cmap=cmap)
+ return hv_plot
+
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):
@@ -674,7 +727,7 @@ def normalize_components(self):
if self.YrA is not None:
self.YrA = nA_mat * self.YrA
if self.R is not None:
- self.R = nA_mat * self.R
+ self.R = nA_mat * self.YrA
if self.bl is not None:
self.bl = nA * self.bl
if self.c1 is not None:
@@ -834,7 +887,7 @@ def evaluate_components(self, imgs, params, dview=None):
CNN classifier threshold
Returns:
- self: esimates object
+ self: estimates object
self.idx_components: np.array
indices of accepted components
self.idx_components_bad: np.array
@@ -849,19 +902,18 @@ def evaluate_components(self, imgs, params, dview=None):
dims = imgs.shape[1:]
opts = params.get_group('quality')
idx_components, idx_components_bad, SNR_comp, r_values, cnn_preds = \
- estimate_components_quality_auto(imgs, self.A, self.C, self.b, self.f,
- self.YrA,
- params.get('data', 'fr'),
- params.get('data', 'decay_time'),
- params.get('init', 'gSig'),
- dims, dview=dview,
- min_SNR=opts['min_SNR'],
- r_values_min=opts['rval_thr'],
- use_cnn=opts['use_cnn'],
- thresh_cnn_min=opts['min_cnn_thr'],
- thresh_cnn_lowest=opts['cnn_lowest'],
- r_values_lowest=opts['rval_lowest'],
- min_SNR_reject=opts['SNR_lowest'])
+ estimate_components_quality_auto(imgs, self.A, self.C, self.b, self.f, self.YrA,
+ params.get('data', 'fr'),
+ params.get('data', 'decay_time'),
+ params.get('init', 'gSig'),
+ dims, dview=dview,
+ min_SNR=opts['min_SNR'],
+ r_values_min=opts['rval_thr'],
+ use_cnn=opts['use_cnn'],
+ thresh_cnn_min=opts['min_cnn_thr'],
+ thresh_cnn_lowest=opts['cnn_lowest'],
+ r_values_lowest=opts['rval_lowest'],
+ min_SNR_reject=opts['SNR_lowest'])
self.idx_components = idx_components.astype(int)
self.idx_components_bad = idx_components_bad.astype(int)
self.SNR_comp = SNR_comp
@@ -962,6 +1014,80 @@ def filter_components(self, imgs, params, new_dict={}, dview=None, select_mode='
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)
+ traces can be performed by setting dff_flag=True. In this case the
+ results of the deconvolution are stored in F_dff_dec and S_dff
+
+ Args:
+ params: params object
+ Parameters of the algorithm
+ dff_flag: bool (True)
+ Flag for deconvolving the DF/F traces
+
+ Returns:
+ self: estimates object
+ '''
+
+ F = self.C + self.YrA
+ args = dict()
+ args['p'] = params.get('preprocess', 'p')
+ args['method_deconvolution'] = params.get('temporal', 'method_deconvolution')
+ args['bas_nonneg'] = params.get('temporal', 'bas_nonneg')
+ args['noise_method'] = params.get('temporal', 'noise_method')
+ args['s_min'] = params.get('temporal', 's_min')
+ args['optimize_g'] = params.get('temporal', 'optimize_g')
+ args['noise_range'] = params.get('temporal', 'noise_range')
+ args['fudge_factor'] = params.get('temporal', 'fudge_factor')
+
+ args_in = [(F[jj], None, jj, None, None, None, None,
+ args) for jj in range(F.shape[0])]
+
+ if 'multiprocessing' in str(type(dview)):
+ results = dview.map_async(
+ constrained_foopsi_parallel, args_in).get(4294967)
+ elif dview is not None:
+ results = dview.map_sync(constrained_foopsi_parallel, args_in)
+ else:
+ results = list(map(constrained_foopsi_parallel, args_in))
+
+ results = list(zip(*results))
+
+ order = list(results[7])
+ self.C = np.stack([results[0][i] for i in order])
+ self.S = np.stack([results[1][i] for i in order])
+ self.bl = [results[3][i] for i in order]
+ self.c1 = [results[4][i] for i in order]
+ self.g = [results[6][i] for i in order]
+ self.neurons_sn = [results[5][i] for i in order]
+ self.lam = [results[8][i] for i in order]
+ self.YrA = F - self.C
+
+ if dff_flag:
+ if self.F_dff is None:
+ logging.warning('The F_dff field is empty. Run the method' +
+ ' estimates.detrend_df_f before attempting' +
+ ' to deconvolve.')
+ else:
+ args_in = [(self.F_dff[jj], None, jj, 0, 0, self.g[jj], None,
+ args) for jj in range(F.shape[0])]
+
+ if 'multiprocessing' in str(type(dview)):
+ results = dview.map_async(
+ constrained_foopsi_parallel, args_in).get(4294967)
+ elif dview is not None:
+ results = dview.map_sync(constrained_foopsi_parallel,
+ args_in)
+ else:
+ results = list(map(constrained_foopsi_parallel, args_in))
+
+ results = list(zip(*results))
+ order = list(results[7])
+ self.F_dff_dec = np.stack([results[0][i] for i in order])
+ self.S_dff = np.stack([results[1][i] for i in order])
+
+
def manual_merge(self, components, params):
''' merge a given list of components. The indices
of components are not pythonic, i.e., they start from 1. Moreover,
@@ -998,6 +1124,7 @@ def manual_merge(self, components, params):
# we initialize the values
A_merged = scipy.sparse.lil_matrix((d, nbmrg))
C_merged = np.zeros((nbmrg, T))
+ R_merged = np.zeros((nbmrg, T))
S_merged = np.zeros((nbmrg, T))
bl_merged = np.zeros((nbmrg, 1))
c1_merged = np.zeros((nbmrg, 1))
@@ -1019,11 +1146,12 @@ def manual_merge(self, components, params):
g_idx = [merged_ROI[indx]]
fast_merge = True
bm, cm, computedA, computedC, gm, \
- sm, ss = merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge,
- None, g_idx, indx, params.temporal)
+ sm, ss, yra = merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge,
+ None, g_idx, indx, params.temporal)
- A_merged[:, i] = computedA
+ A_merged[:, i] = computedA[:, np.newaxis]
C_merged[i, :] = computedC
+ R_merged[i, :] = yra
S_merged[i, :] = ss[:T]
bl_merged[i] = bm
c1_merged[i] = cm
@@ -1035,6 +1163,7 @@ def manual_merge(self, components, params):
if np.any(empty):
A_merged = A_merged[:, ~empty]
C_merged = C_merged[~empty]
+ R_merged = R_merged[~empty]
S_merged = S_merged[~empty]
bl_merged = bl_merged[~empty]
c1_merged = c1_merged[~empty]
@@ -1062,6 +1191,9 @@ def manual_merge(self, components, params):
A_merged.tocsc()))
self.C = np.vstack((self.C[good_neurons, :], C_merged))
# we continue for the variables
+ if self.YrA is not None:
+ self.YrA = np.vstack((self.YrA[good_neurons, :], R_merged))
+ self.R = self.YrA
if self.S is not None:
self.S = np.vstack((self.S[good_neurons, :], S_merged))
if self.bl is not None:
@@ -1080,7 +1212,7 @@ def manual_merge(self, components, params):
self.coordinates = caiman.utils.visualization.get_contours(self.A,\
self.dims, thr_method='max', thr='0.2')
- def threshold_spatial_components(self, maxthr=0.25, dview=None, **kwargs):
+ def threshold_spatial_components(self, maxthr=0.25, dview=None):
''' threshold spatial components. See parameters of
spatial.threshold_components
@@ -1095,8 +1227,10 @@ def threshold_spatial_components(self, maxthr=0.25, dview=None, **kwargs):
'''
if self.A_thr is None:
- A_thr = threshold_components(self.A, self.dims, maxthr=maxthr, dview=dview, medw=None, thr_method='max', nrgthr=0.99,
- extract_cc=True, se=None, ss=None, **kwargs)
+ A_thr = threshold_components(self.A, self.dims, maxthr=maxthr, dview=dview,
+ medw=None, thr_method='max', nrgthr=0.99,
+ extract_cc=True, se=None, ss=None)
+
self.A_thr = A_thr
else:
print('A_thr already computed. If you want to recompute set self.A_thr to None')
@@ -1195,89 +1329,186 @@ def save_NWB(self,
imaging_series_name=None,
sess_desc='CaImAn Results',
exp_desc=None,
- imaging_rate=30,
- location='somewhere in the brain',
- orig_file_format='tiff'):
- """save object in hdf5 file format
+ identifier=None,
+ imaging_rate=30.,
+ starting_time=0.,
+ session_start_time=None,
+ excitation_lambda=488.0,
+ imaging_plane_description='some imaging plane description',
+ emission_lambda=520.0,
+ indicator='OGB-1',
+ location='brain',
+ raw_data_file=None):
+ """writes NWB file
Args:
filename: str
- path to the hdf5 file containing the saved object
+
+ imaging_plane_name: str, optional
+
+ imaging_series_name: str, optional
+
+ sess_desc: str, optional
+
+ exp_desc: str, optional
+
+ identifier: str, optional
+
+ imaging_rate: float, optional
+ default: 30 (Hz)
+
+ starting_time: float, optional
+ default: 0.0 (seconds)
+
+ location: str, optional
+
+ session_start_time: datetime.datetime, optional
+ Only required for new files
+
+ excitation_lambda: float
+
+ imaging_plane_description: str
+
+ emission_lambda: float
+
+ indicator: str
+
+ location: str
"""
- from pynwb import NWBHDF5IO
- from pynwb.ophys import ImageSegmentation, Fluorescence, MotionCorrection
+
+ from pynwb import NWBHDF5IO, TimeSeries, NWBFile
+ from pynwb.base import Images
+ from pynwb.image import GrayscaleImage
+ from pynwb.ophys import ImageSegmentation, Fluorescence, OpticalChannel, ImageSeries
+ from pynwb.device import Device
import os
+
+ if identifier is None:
+ import uuid
+ identifier = uuid.uuid1().hex
+
if '.nwb' != os.path.splitext(filename)[-1].lower():
raise Exception("Wrong filename")
- if not os.path.isfile(filename): # if the file doesn't exist create new and add the orginal data path
- raise Exception('filename should be an existing NWB file.\
- Consider using the cnmf.movie.save method to create one.')
-
- else: # if the file already exist in the .nwb format then just add the results to it
- logging.info('Saving the results in the NWB file...')
- with NWBHDF5IO(filename, 'r+') as io:
- nwbfile = io.read()
- # Add processing results
- mod = nwbfile.create_processing_module('ophys', 'contains caiman estimates for the main imagin plane')
- img_seg = ImageSegmentation()
- mod.add_data_interface(img_seg)
- fl = Fluorescence()
- mod.add_data_interface(fl)
- # mot_crct = MotionCorrection()
- # mod.add_data_interface(mot_crct)
-
- # Add the ROI-related stuff
- if imaging_plane_name is None:
- imaging_plane_name = [imp for imp in nwbfile.imaging_planes.keys()]
- if len(imaging_plane_name)>1:
- raise Exception('There is more than one imaging plane in the file, you need to specify the name via '
- 'the "imaging_plane_name" parameter')
- else:
- imaging_plane_name = imaging_plane_name[0]
+ if not os.path.isfile(filename): # if the file doesn't exist create new and add the original data path
+ print('filename does not exist. Creating new NWB file with only estimates output')
+
+ nwbfile = NWBFile(sess_desc, identifier, session_start_time, experiment_description=exp_desc)
+ device = Device('imaging_device')
+ nwbfile.add_device(device)
+ optical_channel = OpticalChannel('OpticalChannel',
+ 'main optical channel',
+ emission_lambda=emission_lambda)
+ nwbfile.create_imaging_plane(name='ImagingPlane',
+ optical_channel=optical_channel,
+ description=imaging_plane_description,
+ device=device,
+ excitation_lambda=excitation_lambda,
+ imaging_rate=imaging_rate,
+ indicator=indicator,
+ location=location)
+ if raw_data_file:
+ nwbfile.add_acquisition(ImageSeries(name='TwoPhotonSeries',
+ external_file=[raw_data_file],
+ format='external',
+ rate=imaging_rate,
+ starting_frame=[0]))
+ with NWBHDF5IO(filename, 'w') as io:
+ io.write(nwbfile)
- if imaging_series_name is None:
- imaging_series_name = [imp for imp in nwbfile.acquisition.keys()]
- if len(imaging_series_name)>1:
- raise Exception('There is more than one imaging plane in the file, you need to specify the name via '
- 'the "imaging_series_name" parameter')
- else:
- imaging_series_name = imaging_series_name[0]
+ time.sleep(4) # ensure the file is fully closed before opening again in append mode
+ logging.info('Saving the results in the NWB file...')
+
+ with NWBHDF5IO(filename, 'r+') as io:
+ nwbfile = io.read()
+ # Add processing results
+
+ # Create the module as 'ophys' unless it is taken and append 'ophysX' instead
+ ophysmodules = [s[5:] for s in list(nwbfile.modules) if s.startswith('ophys')]
+ if any('' in s for s in ophysmodules):
+ if any([s for s in ophysmodules if s.isdigit()]):
+ nummodules = max([int(s) for s in ophysmodules if s.isdigit()])+1
+ print('ophys module previously created, writing to ophys'+str(nummodules)+' instead')
+ mod = nwbfile.create_processing_module('ophys'+str(nummodules), 'contains caiman estimates for '
+ 'the main imaging plane')
+ else:
+ print('ophys module previously created, writing to ophys1 instead')
+ mod = nwbfile.create_processing_module('ophys1', 'contains caiman estimates for the main '
+ 'imaging plane')
+ else:
+ mod = nwbfile.create_processing_module('ophys', 'contains caiman estimates for the main imaging plane')
+ img_seg = ImageSegmentation()
+ mod.add_data_interface(img_seg)
+ fl = Fluorescence()
+ mod.add_data_interface(fl)
+# mot_crct = MotionCorrection()
+# mod.add_data_interface(mot_crct)
+ # Add the ROI-related stuff
+ if imaging_plane_name is not None:
imaging_plane = nwbfile.imaging_planes[imaging_plane_name]
- image_series = nwbfile.acquisition[imaging_series_name]
-
- ps = img_seg.create_plane_segmentation('CNMF_ROIs',
- imaging_plane, 'planeseg', image_series)
-
- # Add ROIs
- for roi in self.A.T: # Neurons
- ps.add_roi(image_mask=roi.T.toarray().reshape(self.dims))
- for bg in self.b.T: # Backgrounds
- ps.add_roi(image_mask=bg.reshape(self.dims))
- # Add Traces
- n_rois = self.A.shape[-1]
- n_bg = len(self.f)
- rt_region_roi = ps.create_roi_table_region('ROIs',
- region=list(range(n_rois)))
-
- rt_region_bg = ps.create_roi_table_region('Background',
- region=list(range(n_rois,n_rois+n_bg)))
-
- timestamps = list(range(self.f.shape[1]))
+ else:
+ if len(nwbfile.imaging_planes) == 1:
+ imaging_plane = list(nwbfile.imaging_planes.values())[0]
+ else:
+ raise Exception('There is more than one imaging plane in the file, you need to specify the name'
+ ' via the "imaging_plane_name" parameter')
- # Neurons
- rrs1 = fl.create_roi_response_series('RoiResponseSeries', self.C.T, 'lumens', rt_region_roi, timestamps=timestamps)
- # Background
- rrs2 = fl.create_roi_response_series('Background_Fluorescence_Response', self.f.T, 'lumens', rt_region_bg, timestamps=timestamps)
+ if imaging_series_name is not None:
+ image_series = nwbfile.acquisition[imaging_series_name]
+ else:
+ if not len(nwbfile.acquisition):
+ image_series = None
+ elif len(nwbfile.acquisition) == 1:
+ image_series = list(nwbfile.acquisition.values())[0]
+ else:
+ raise Exception('There is more than one imaging plane in the file, you need to specify the name'
+ ' via the "imaging_series_name" parameter')
+
+ ps = img_seg.create_plane_segmentation('CNMF_ROIs', imaging_plane, 'PlaneSegmentation', image_series)
+
+ ps.add_column('r', 'description of r values')
+ ps.add_column('snr', 'signal to noise ratio')
+ ps.add_column('cnn', 'description of CNN')
+ ps.add_column('keep', 'in idx_components')
+
+ # Add ROIs
+ for i, (roi, snr, r, cnn) in enumerate(zip(self.A.T, self.SNR_comp, self.r_values, self.cnn_preds)):
+ ps.add_roi(image_mask=roi.T.toarray().reshape(self.dims), r=r, snr=snr, cnn=cnn,
+ keep=i in self.idx_components)
+ for bg in self.b.T: # Backgrounds
+ ps.add_roi(image_mask=bg.reshape(self.dims), r=np.nan, snr=np.nan, cnn=np.nan, keep=False)
+ # Add Traces
+ n_rois = self.A.shape[-1]
+ n_bg = len(self.f)
+ rt_region_roi = ps.create_roi_table_region(
+ 'ROIs', region=list(range(n_rois)))
+
+ rt_region_bg = ps.create_roi_table_region(
+ 'Background', region=list(range(n_rois, n_rois+n_bg)))
+
+ timestamps = np.arange(self.f.shape[1]) / imaging_rate + starting_time
+
+ # Neurons
+ fl.create_roi_response_series('RoiResponseSeries', self.C.T, 'lumens', rt_region_roi, timestamps=timestamps)
+ # Background
+ fl.create_roi_response_series('Background_Fluorescence_Response', self.f.T, 'lumens', rt_region_bg,
+ timestamps=timestamps)
+
+ mod.add(TimeSeries(name='residuals', description='residuals', data=self.YrA.T, timestamps=timestamps,
+ unit='NA'))
+ if hasattr(self, 'Cn'):
+ images = Images('summary_images')
+ images.add_image(GrayscaleImage(name='local_correlations', data=self.Cn))
# Add MotionCorreciton
# create_corrected_image_stack(corrected, original, xy_translation, name='CorrectedImageStack')
io.write(nwbfile)
-def compare_components(estimate_gt, estimate_cmp, Cn=None, thresh_cost=.8, min_dist=10, print_assignment=False, labels=['GT', 'CMP'], plot_results=False):
+def compare_components(estimate_gt, estimate_cmp, Cn=None, thresh_cost=.8, min_dist=10, print_assignment=False,
+ labels=['GT', 'CMP'], plot_results=False):
if estimate_gt.A_thr is None:
raise Exception(
'You need to compute thresolded components for first argument before calling remove_duplicates: use the threshold_components method')
diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py
index d69902e90..fa0e46dc4 100644
--- a/caiman/source_extraction/cnmf/initialization.py
+++ b/caiman/source_extraction/cnmf/initialization.py
@@ -38,6 +38,7 @@
import caiman
from .deconvolution import constrained_foopsi
+#from .utilities import fast_graph_Laplacian_patches
from .pre_processing import get_noise_fft, get_noise_welch
from .spatial import circular_constraint, connectivity_constraint
@@ -140,7 +141,9 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter
max_iter_snmf=500, alpha_snmf=10e2, sigma_smooth_snmf=(.5, .5, .5),
perc_baseline_snmf=20, options_local_NMF=None, rolling_sum=False,
rolling_length=100, sn=None, options_total=None, min_corr=0.8, min_pnr=10,
- ring_size_factor=1.5, center_psf=False, ssub_B=2, init_iter=2, remove_baseline = True):
+ 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):
"""
Initalize components
@@ -319,6 +322,14 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter
Ain, Cin, _, b_in, f_in = sparseNMF(
Y_ds, nr=K, nb=nb, max_iter_snmf=max_iter_snmf, alpha=alpha_snmf,
sigma_smooth=sigma_smooth_snmf, remove_baseline=remove_baseline, perc_baseline=perc_baseline_snmf)
+
+ elif method == 'graph_nmf':
+ Ain, Cin, _, b_in, f_in = graphNMF(
+ Y_ds, nr=K, nb=nb, max_iter_snmf=max_iter_snmf, lambda_gnmf=lambda_gnmf,
+ sigma_smooth=sigma_smooth_snmf, remove_baseline=remove_baseline,
+ perc_baseline=perc_baseline_snmf, SC_kernel=SC_kernel,
+ 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':
Ain, Cin, _, b_in, f_in = ICA_PCA(
@@ -475,7 +486,6 @@ def ICA_PCA(Y_ds, nr, sigma_smooth=(.5, .5, .5), truncate=2, fun='logcosh',
model = NMF(n_components=nb, init='random', random_state=0)
-
b_in = model.fit_transform(np.maximum(m1, 0)).astype(np.float32)
f_in = model.components_.astype(np.float32)
@@ -489,6 +499,12 @@ def sparseNMF(Y_ds, nr, max_iter_snmf=500, alpha=10e2, sigma_smooth=(.5, .5, .5)
Initialization using sparse NMF
Args:
+ Y_ds: nd.array or movie (T, x, y [,z])
+ data
+
+ nr: int
+ number of components
+
max_iter_snm: int
number of iterations
@@ -506,8 +522,8 @@ def sparseNMF(Y_ds, nr, max_iter_snmf=500, alpha=10e2, sigma_smooth=(.5, .5, .5)
Returns:
A: np.array
- 2d array of size (# of pixels) x nr with the spatial components. Each column is
- ordered columnwise (matlab format, order='F')
+ 2d array of size (# of pixels) x nr with the spatial components.
+ Each column is ordered columnwise (matlab format, order='F')
C: np.array
2d array of size nr X T with the temporal components
@@ -517,7 +533,8 @@ def sparseNMF(Y_ds, nr, max_iter_snmf=500, alpha=10e2, sigma_smooth=(.5, .5, .5)
"""
m = scipy.ndimage.gaussian_filter(np.transpose(
- Y_ds, [2, 0, 1]), sigma=sigma_smooth, mode='nearest', truncate=truncate)
+ Y_ds, np.roll(np.arange(Y_ds.ndim), 1)), sigma=sigma_smooth,
+ mode='nearest', truncate=truncate)
if remove_baseline:
logging.info('REMOVING BASELINE')
bl = np.percentile(m, perc_baseline, axis=0)
@@ -527,11 +544,11 @@ def sparseNMF(Y_ds, nr, max_iter_snmf=500, alpha=10e2, sigma_smooth=(.5, .5, .5)
bl = np.zeros(m.shape[1:])
m1 = m
- T, d1, d2 = np.shape(m1)
- d = d1 * d2
+ T, dims = m1.shape[0], m1.shape[1:]
+ d = np.prod(dims)
yr = np.reshape(m1, [T, d], order='F')
mdl = NMF(n_components=nr, verbose=False, init='nndsvd', tol=1e-10,
- max_iter=max_iter_snmf, shuffle=False, alpha=alpha, l1_ratio=1)
+ max_iter=max_iter_snmf, shuffle=False, alpha=alpha, l1_ratio=1)
C = mdl.fit_transform(yr).T
A = mdl.components_.T
A_in = A
@@ -542,7 +559,59 @@ def sparseNMF(Y_ds, nr, max_iter_snmf=500, alpha=10e2, sigma_smooth=(.5, .5, .5)
random_state=0, max_iter=max_iter_snmf)
b_in = model.fit_transform(np.maximum(m1, 0)).astype(np.float32)
f_in = model.components_.astype(np.float32)
- center = caiman.base.rois.com(A_in, d1, d2)
+ center = caiman.base.rois.com(A_in, *dims)
+
+ return A_in, C_in, center, b_in, f_in
+
+def graphNMF(Y_ds, nr, max_iter_snmf=500, lambda_gnmf=1,
+ sigma_smooth=(.5, .5, .5), remove_baseline=True,
+ perc_baseline=20, nb=1, truncate=2, tol=1e-3, SC_kernel='heat',
+ SC_normalize=True, SC_thr=0, SC_sigma=1, SC_use_NN=False,
+ SC_nnn=20):
+
+ m = scipy.ndimage.gaussian_filter(np.transpose(
+ Y_ds, np.roll(np.arange(Y_ds.ndim), 1)), sigma=sigma_smooth,
+ mode='nearest', truncate=truncate)
+ if remove_baseline:
+ logging.info('REMOVING BASELINE')
+ bl = np.percentile(m, perc_baseline, axis=0)
+ m1 = np.maximum(0, m - bl)
+ else:
+ logging.info('NOT REMOVING BASELINE')
+ bl = np.zeros(m.shape[1:])
+ m1 = m
+
+ T, dims = m1.shape[0], m1.shape[1:]
+ d = np.prod(dims)
+ yr = np.reshape(m1, [T, d], order='F')
+ mdl = NMF(n_components=nr, verbose=False, init='nndsvd', tol=1e-10,
+ max_iter=5)
+ C = mdl.fit_transform(yr).T
+ A = mdl.components_.T
+ W = caiman.source_extraction.cnmf.utilities.fast_graph_Laplacian_patches(
+ [np.reshape(m, [T, d], order='F').T, [], 'heat', SC_sigma, SC_thr,
+ SC_nnn, SC_normalize, SC_use_NN])
+ D = scipy.sparse.spdiags(W.sum(0), 0, W.shape[0], W.shape[0])
+ for it in range(max_iter_snmf):
+ C_ = C.copy()
+ A_ = A.copy()
+ C = C*(yr.dot(A)/(C.T.dot(A.T.dot(A))+np.finfo(C.dtype).eps)).T
+ A = A*(yr.T.dot(C.T) + lambda_gnmf*(W.dot(A)))/(A.dot(C.dot(C.T)) + lambda_gnmf*D.dot(A) + np.finfo(C.dtype).eps)
+ nA = np.sqrt((A**2).sum(0))
+ A /= nA
+ C *= nA[:, np.newaxis]
+ if (np.linalg.norm(C - C_)/np.linalg.norm(C_) < tol) & (np.linalg.norm(A - A_)/np.linalg.norm(A_) < tol):
+ logging.info('Graph NMF converged after {} iterations'.format(it+1))
+ break
+ A_in = A
+ C_in = C
+
+ m1 = yr.T - A_in.dot(C_in) + np.maximum(0, bl.flatten(order='F'))[:, np.newaxis]
+ model = NMF(n_components=nb, init='random',
+ random_state=0, max_iter=max_iter_snmf)
+ b_in = model.fit_transform(np.maximum(m1, 0)).astype(np.float32)
+ f_in = model.components_.astype(np.float32)
+ center = caiman.base.rois.com(A_in, *dims)
return A_in, C_in, center, b_in, f_in
@@ -999,7 +1068,7 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p
# 1st iteration on decimated data
logging.info('Merging components')
A, C = caiman.source_extraction.cnmf.merging.merge_components(
- B, A, [], C, [], C, [], o, options['spatial_params'],
+ B, A, [], C, None, [], C, [], o, options['spatial_params'],
dview=None, thr=options['merging']['merge_thr'], mx=np.Inf, fast_merge=True)[:2]
A = A.astype(np.float32)
C = C.astype(np.float32)
@@ -1048,7 +1117,7 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p
logging.info('Merging components')
A, C = caiman.source_extraction.cnmf.merging.merge_components(
- B, A, [], C, [], C, [], o, options['spatial_params'],
+ B, A, [], C, None, [], C, [], o, options['spatial_params'],
dview=None, thr=options['merging']['merge_thr'], mx=np.Inf, fast_merge=True)[:2]
A = A.astype(np.float32)
C = C.astype(np.float32)
diff --git a/caiman/source_extraction/cnmf/merging.py b/caiman/source_extraction/cnmf/merging.py
index 852c274de..20357fcbd 100644
--- a/caiman/source_extraction/cnmf/merging.py
+++ b/caiman/source_extraction/cnmf/merging.py
@@ -15,7 +15,7 @@
import logging
from past.utils import old_div
import scipy
-from scipy.sparse import coo_matrix, csgraph, csc_matrix, lil_matrix
+from scipy.sparse import csgraph, csc_matrix, lil_matrix, csr_matrix
from .spatial import update_spatial_components, threshold_components
from .temporal import update_temporal_components
@@ -24,51 +24,61 @@
-def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params, dview=None, thr=0.85, fast_merge=True, mx=1000, bl=None, c1=None, sn=None, g=None):
+def merge_components(Y, A, b, C, R, f, S, sn_pix, temporal_params,
+ spatial_params, dview=None, thr=0.85, fast_merge=True,
+ mx=1000, bl=None, c1=None, sn=None, g=None,
+ merge_parallel=False, max_merge_area=None):
+
""" Merging of spatially overlapping components that have highly correlated temporal activity
The correlation threshold for merging overlapping components is user specified in thr
Args:
Y: np.ndarray
- residual movie after subtracting all found components (Y_res = Y - A*C - b*f) (d x T)
+ residual movie after subtracting all found components
+ (Y_res = Y - A*C - b*f) (d x T)
A: sparse matrix
matrix of spatial components (d x K)
b: np.ndarray
spatial background (vector of length d)
-
+
C: np.ndarray
matrix of temporal components (K x T)
-
+
+ R: np.ndarray
+ array of residuals (K x T)
+
f: np.ndarray
temporal background (vector of length T)
-
+
S: np.ndarray
matrix of deconvolved activity (spikes) (K x T)
-
+
sn_pix: ndarray
noise standard deviation for each pixel
-
+
temporal_params: dictionary
- all the parameters that can be passed to the update_temporal_components function
-
+ all the parameters that can be passed to the
+ update_temporal_components function
+
spatial_params: dictionary
- all the parameters that can be passed to the update_spatial_components function
-
+ all the parameters that can be passed to the
+ update_spatial_components function
+
thr: scalar between 0 and 1
correlation threshold for merging (default 0.85)
-
+
mx: int
maximum number of merging operations (default 50)
-
+
sn_pix: nd.array
noise level for each pixel (vector of length d)
-
+
fast_merge: bool
if true perform rank 1 merging, otherwise takes best neuron
-
+
bl:
baseline for fluorescence trace for each row in C
c1:
@@ -78,40 +88,50 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params,
sn:
noise level for each row in C
+ merge_parallel: bool
+ perform merging in parallel
+
+ max_merge_area: int
+ maximum area (in pixels) of merged components,
+ used to determine whether to merge
+
Returns:
A: sparse matrix
matrix of merged spatial components (d x K)
-
+
C: np.ndarray
matrix of merged temporal components (K x T)
-
+
nr: int
number of components after merging
-
+
merged_ROIs: list
index of components that have been merged
-
+
S: np.ndarray
matrix of merged deconvolved activity (spikes) (K x T)
-
+
bl: float
baseline for fluorescence trace
-
+
c1: float
initial concentration
-
+
g: float
discrete time constant
-
+
sn: float
noise level
+ R: np.ndarray
+ residuals
Raises:
Exception "The number of elements of bl\c1\g\sn must match the number of components"
"""
#tests and initialization
nr = A.shape[1]
+ A = csc_matrix(A)
if bl is not None and len(bl) != nr:
raise Exception(
"The number of elements of bl must match the number of components")
@@ -124,6 +144,8 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params,
if g is not None and len(g) != nr:
raise Exception(
"The number of elements of g must match the number of components")
+ if R is None:
+ R = np.zeros_like(C)
[d, t] = np.shape(Y)
@@ -172,114 +194,157 @@ def merge_components(Y, A, b, C, f, S, sn_pix, temporal_params, spatial_params,
nbmrg = min((np.size(ind), mx)) # number of merging operations
- # we initialize the values
- A_merged = lil_matrix((d, nbmrg))
- C_merged = np.zeros((nbmrg, t))
- S_merged = np.zeros((nbmrg, t))
- bl_merged = np.zeros((nbmrg, 1))
- c1_merged = np.zeros((nbmrg, 1))
- sn_merged = np.zeros((nbmrg, 1))
- g_merged = np.zeros((nbmrg, p))
- merged_ROIs = []
-
- for i in range(nbmrg):
- merged_ROI = np.where(list_conxcomp[:, ind[i]])[0]
- logging.info('Merging components {}'.format(merged_ROI))
- merged_ROIs.append(merged_ROI)
-
- Acsc = A.tocsc()[:, merged_ROI]
- Ctmp = np.array(C)[merged_ROI, :]
-
-
- # # we l2 the traces to have normalization values
- # C_to_norm = np.sqrt([computedC.dot(computedC)
- # for computedC in C[merged_ROI]])
-# fast_merge = False
-
- # from here we are computing initial values for C and A
-
-
- # this is a big normalization value that for every one of the merged neuron
- C_to_norm = np.sqrt(np.ravel(Acsc.power(2).sum(
- axis=0)) * np.sum(Ctmp ** 2, axis=1))
- indx = np.argmax(C_to_norm)
- g_idx = [merged_ROI[indx]]
-
- bm, cm, computedA, computedC, gm, sm, ss = merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx,
- indx, temporal_params)
-
- A_merged[:, i] = computedA
- C_merged[i, :] = computedC
- S_merged[i, :] = ss[:t]
- bl_merged[i] = bm
- c1_merged[i] = cm
- sn_merged[i] = sm
- g_merged[i, :] = gm
+ if merge_parallel:
+ merged_ROIs = [np.where(list_conxcomp[:, ind[i]])[0] for i in range(nbmrg)]
+ Acsc_mats = [csc_matrix(A[:, merged_ROI]) for merged_ROI in merged_ROIs]
+ Ctmp_mats = [C[merged_ROI] + R[merged_ROI] for merged_ROI in merged_ROIs]
+ C_to_norms = [np.sqrt(np.ravel(Acsc.power(2).sum(
+ axis=0)) * np.sum(Ctmp ** 2, axis=1)) for (Acsc, Ctmp) in zip(Acsc_mats, Ctmp_mats)]
+ indxs = [np.argmax(C_to_norm) for C_to_norm in C_to_norms]
+ g_idxs = [merged_ROI[indx] for (merged_ROI, indx) in zip(merged_ROIs, indxs)]
+ fms = [fast_merge]*nbmrg
+ tps = [temporal_params]*nbmrg
+ gs = [g]*nbmrg
+
+ if dview is None:
+ merge_res = list(map(merge_iter, zip(Acsc_mats, C_to_norms, Ctmp_mats, fms, gs, g_idxs, indxs, tps)))
+ elif 'multiprocessin' in str(type(dview)):
+ merge_res = list(dview.map(merge_iter, zip(Acsc_mats, C_to_norms, Ctmp_mats, fms, gs, g_idxs, indxs, tps)))
+ else:
+ merge_res = list(dview.map_sync(merge_iter, zip(Acsc_mats, C_to_norms, Ctmp_mats, fms, gs, g_idxs, indxs, tps)))
+ dview.results.clear()
+ #merge_res = list(dview.map(merge_iter, zip(Acsc_mats, C_to_norms, Ctmp_mats, fms, gs, g_idxs, indxs, tps)))
+ bl_merged = np.array([res[0] for res in merge_res])
+ c1_merged = np.array([res[1] for res in merge_res])
+ A_merged = scipy.sparse.vstack([csc_matrix(res[2]) for res in merge_res]).T
+ C_merged = np.vstack([res[3] for res in merge_res])
+ g_merged = np.vstack([res[4] for res in merge_res])
+ sn_merged = np.array([res[5] for res in merge_res])
+ S_merged = np.vstack([res[6] for res in merge_res])
+ R_merged = np.vstack([res[7] for res in merge_res])
+ else:
+ # we initialize the values
+ A_merged = lil_matrix((d, nbmrg))
+ C_merged = np.zeros((nbmrg, t))
+ R_merged = np.zeros((nbmrg, t))
+ S_merged = np.zeros((nbmrg, t))
+ bl_merged = np.zeros((nbmrg, 1))
+ c1_merged = np.zeros((nbmrg, 1))
+ sn_merged = np.zeros((nbmrg, 1))
+ g_merged = np.zeros((nbmrg, p))
+ merged_ROIs = []
+ for i in range(nbmrg):
+ merged_ROI = np.where(list_conxcomp[:, ind[i]])[0]
+ logging.info('Merging components {}'.format(merged_ROI))
+ merged_ROIs.append(merged_ROI)
+ Acsc = A.tocsc()[:, merged_ROI]
+ Ctmp = np.array(C)[merged_ROI, :] + np.array(R)[merged_ROI, :]
+ C_to_norm = np.sqrt(np.ravel(Acsc.power(2).sum(
+ axis=0)) * np.sum(Ctmp ** 2, axis=1))
+ indx = np.argmax(C_to_norm)
+ g_idx = [merged_ROI[indx]]
+ bm, cm, computedA, computedC, gm, sm, ss, yra = merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx,
+ indx, temporal_params)
+
+ A_merged[:, i] = csr_matrix(computedA).T
+ C_merged[i, :] = computedC
+ R_merged[i, :] = yra
+ S_merged[i, :] = ss[:t]
+ bl_merged[i] = bm
+ c1_merged[i] = cm
+ sn_merged[i] = sm
+ g_merged[i, :] = gm
empty = np.ravel((C_merged.sum(1) == 0) + (A_merged.sum(0) == 0))
if np.any(empty):
A_merged = A_merged[:, ~empty]
C_merged = C_merged[~empty]
+ R_merged = R_merged[~empty]
S_merged = S_merged[~empty]
bl_merged = bl_merged[~empty]
c1_merged = c1_merged[~empty]
sn_merged = sn_merged[~empty]
g_merged = g_merged[~empty]
- # we want to remove merged neuron from the initial part and replace them with merged ones
- neur_id = np.unique(np.hstack(merged_ROIs))
- good_neurons = np.setdiff1d(list(range(nr)), neur_id)
- A = scipy.sparse.hstack((A.tocsc()[:, good_neurons], A_merged.tocsc()))
- C = np.vstack((C[good_neurons, :], C_merged))
- # we continue for the variables
- if S is not None:
- S = np.vstack((S[good_neurons, :], S_merged))
- if bl is not None:
- bl = np.hstack((bl[good_neurons], np.array(bl_merged).flatten()))
- if c1 is not None:
- c1 = np.hstack((c1[good_neurons], np.array(c1_merged).flatten()))
- if sn is not None:
- sn = np.hstack((sn[good_neurons], np.array(sn_merged).flatten()))
- if g is not None:
- g = np.vstack((np.vstack(g)[good_neurons], g_merged))
- nr = nr - len(neur_id) + len(C_merged)
+ if len(merged_ROIs) > 0:
+ # we want to remove merged neuron from the initial part and replace them with merged ones
+ neur_id = np.unique(np.hstack(merged_ROIs))
+ good_neurons = np.setdiff1d(list(range(nr)), neur_id)
+ A = scipy.sparse.hstack((A.tocsc()[:, good_neurons], A_merged.tocsc()))
+ C = np.vstack((C[good_neurons, :], C_merged))
+ # we continue for the variables
+ if S is not None:
+ S = np.vstack((S[good_neurons, :], S_merged))
+ if R is not None:
+ R = np.vstack((R[good_neurons, :], R_merged))
+ if bl is not None:
+ bl = np.hstack((bl[good_neurons], np.array(bl_merged).flatten()))
+ if c1 is not None:
+ c1 = np.hstack((c1[good_neurons], np.array(c1_merged).flatten()))
+ if sn is not None:
+ sn = np.hstack((sn[good_neurons], np.array(sn_merged).flatten()))
+ if g is not None:
+ g = np.vstack((np.vstack(g)[good_neurons], g_merged))
+# g = np.vstack(g)[good_neurons]
+# if g.shape[1] == 0:
+# g = np.zeros((len(good_neurons), 1))
+# g = np.vstack((g, g_merged))
+
+ nr = nr - len(neur_id) + len(C_merged)
else:
logging.info('No more components merged!')
merged_ROIs = []
empty = []
- return A, C, nr, merged_ROIs, S, bl, c1, sn, g, empty
+ return A, C, nr, merged_ROIs, S, bl, c1, sn, g, empty, R
+def merge_iter(a):
+ Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx, indx, temporal_params = a
+ res = merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx,
+ indx, temporal_params)
+ return res
def merge_iteration(Acsc, C_to_norm, Ctmp, fast_merge, g, g_idx, indx, temporal_params):
if fast_merge:
# we normalize the values of different A's to be able to compare them efficiently. we then sum them
- computedA = Acsc.dot(scipy.sparse.diags(
- C_to_norm, 0, (len(C_to_norm), len(C_to_norm)))).sum(axis=1)
- # we operate a rank one NMF, refining it multiple times (see cnmf demos )
+ computedA = Acsc.dot(C_to_norm)
for _ in range(10):
- computedC = np.maximum(Acsc.T.dot(computedA).T.dot(
- Ctmp) / (computedA.T * computedA), 0)
- if computedC * computedC.T == 0:
+ computedC = np.maximum((Acsc.T.dot(computedA)).dot(Ctmp) /
+ (computedA.T.dot(computedA)), 0)
+ nc = computedC.T.dot(computedC)
+ if nc == 0:
break
- computedA = np.maximum(
- Acsc.dot(Ctmp.dot(computedC.T)) / (computedC * computedC.T), 0)
+ computedA = np.maximum(Acsc.dot(Ctmp.dot(computedC.T)) / nc, 0)
+
+# computedA = Acsc.dot(scipy.sparse.diags(
+# C_to_norm, 0, (len(C_to_norm), len(C_to_norm)))).sum(axis=1)
+#
+# # we operate a rank one NMF, refining it multiple times (see cnmf demos )
+# for _ in range(10):
+# computedC = np.maximum(Acsc.T.dot(computedA).T.dot(
+# Ctmp) / (computedA.T * computedA), 0)
+# if computedC * computedC.T == 0:
+# break
+# computedA = np.maximum(
+# Acsc.dot(Ctmp.dot(computedC.T)) / (computedC * computedC.T), 0)
else:
logging.info('Simple merging ny taking best neuron')
computedC = Ctmp[indx]
computedA = Acsc[:, indx]
# then we de-normalize them using A_to_norm
- A_to_norm = np.sqrt(computedA.T.dot(computedA)[
- 0, 0] / Acsc.power(2).sum(0).max())
+ A_to_norm = np.sqrt(computedA.T.dot(computedA)) #/Acsc.power(2).sum(0).max())
computedA /= A_to_norm
computedC *= A_to_norm
+ # r = (computedA.T.dot(Acsc.dot(Ctmp)))/(computedA.T.dot(computedA)) - computedC
+ r = ((Acsc.T.dot(computedA)).dot(Ctmp))/(computedA.T.dot(computedA)) - computedC
# we then compute the traces ( deconvolution ) to have a clean c and noise in the background
+ c_in = np.array(computedC+r).squeeze()
if g is not None:
- computedC, bm, cm, gm, sm, ss, lam_ = constrained_foopsi(
- np.array(computedC).squeeze(), g=g_idx, **temporal_params)
+ deconvC, bm, cm, gm, sm, ss, lam_ = constrained_foopsi(
+ c_in, g=g_idx, **temporal_params)
else:
- computedC, bm, cm, gm, sm, ss, lam_ = constrained_foopsi(
- np.array(computedC).squeeze(), g=None, **temporal_params)
- return bm, cm, computedA, computedC, gm, sm, ss
+ deconvC, bm, cm, gm, sm, ss, lam_ = constrained_foopsi(
+ c_in, g=None, **temporal_params)
+ return bm, cm, computedA, deconvC, gm, sm, ss, c_in - deconvC
diff --git a/caiman/source_extraction/cnmf/online_cnmf.py b/caiman/source_extraction/cnmf/online_cnmf.py
index 6ac8acd73..fac0d3946 100644
--- a/caiman/source_extraction/cnmf/online_cnmf.py
+++ b/caiman/source_extraction/cnmf/online_cnmf.py
@@ -24,13 +24,14 @@
from past.utils import old_div
from scipy.ndimage import percentile_filter
from scipy.ndimage.filters import gaussian_filter
-from scipy.sparse import coo_matrix, csc_matrix, spdiags
+from scipy.sparse import coo_matrix, csc_matrix, spdiags, hstack
from scipy.stats import norm
from sklearn.decomposition import NMF
from sklearn.preprocessing import normalize
from time import time
from typing import List, Tuple
import tensorflow as tf
+import logging
import caiman
from .cnmf import CNMF
@@ -216,21 +217,18 @@ def _prepare_object(self, Yr, T, new_dims=None, idx_components=None):
# setup per patch classifier
if self.params.get('online', 'path_to_model') is None or self.params.get('online', 'sniper_mode') is False:
- self.loaded_model = None
- self.tf_in = None
- self.tf_out = None
+ loaded_model = None
self.params.set('online', {'sniper_mode': False})
self.tf_in = None
self.tf_out = None
else:
try:
- import keras
- from keras.models import model_from_json
- #logging.debug('Using Keras')
+ from tensorflow.keras.models import model_from_json
+ logging.info('Using Keras')
use_keras = True
except(ModuleNotFoundError):
use_keras = False
- #logging.debug('Using Tensorflow')
+ logging.info('Using Tensorflow')
if use_keras:
path = self.params.get('online', 'path_to_model').split(".")[:-1]
json_path = ".".join(path + ["json"])
@@ -240,9 +238,9 @@ def _prepare_object(self, Yr, T, new_dims=None, idx_components=None):
json_file.close()
loaded_model = model_from_json(loaded_model_json)
loaded_model.load_weights(model_path)
- opt = keras.optimizers.rmsprop(lr=0.0001, decay=1e-6)
- loaded_model.compile(loss=keras.losses.categorical_crossentropy,
- optimizer=opt, metrics=['accuracy'])
+ #opt = tf.keras.optimizers.rmsprop(lr=0.0001, decay=1e-6)
+ #loaded_model.compile(loss=tf.keras.losses.categorical_crossentropy,
+ # optimizer=opt, metrics=['accuracy'])
self.tf_in = None
self.tf_out = None
else:
@@ -831,9 +829,18 @@ def fit_online(self, **kwargs):
if self.estimates.OASISinstances is not None:
self.estimates.bl = [osi.b for osi in self.estimates.OASISinstances]
self.estimates.S = np.stack([osi.s for osi in self.estimates.OASISinstances])
+ self.estimates.S = self.estimates.S[:, t - t // epochs:t]
else:
self.estimates.bl = [0] * self.estimates.C.shape[0]
self.estimates.S = np.zeros_like(self.estimates.C)
+ if self.params.get('online', 'ds_factor') > 1:
+ dims = Y_.shape[1:]
+ self.estimates.A = hstack([coo_matrix(cv2.resize(self.estimates.A[:, i].reshape(self.dims, order='F').toarray(),
+ dims[::-1]).reshape(-1, order='F')[:,None]) for i in range(self.N)], format='csc')
+ self.estimates.b = np.concatenate([cv2.resize(self.estimates.b[:, i].reshape(self.dims, order='F'),
+ dims[::-1]).reshape(-1, order='F')[:,None] for i in range(self.params.get('init', 'nb'))], axis=1)
+ self.params.set('data', {'dims': dims})
+ self.estimates.dims = dims
if self.params.get('online', 'save_online_movie'):
out.release()
if self.params.get('online', 'show_movie'):
diff --git a/caiman/source_extraction/cnmf/params.py b/caiman/source_extraction/cnmf/params.py
index c81bbafcf..08eb65b19 100644
--- a/caiman/source_extraction/cnmf/params.py
+++ b/caiman/source_extraction/cnmf/params.py
@@ -1,10 +1,10 @@
import logging
import os
-
import numpy as np
import scipy
from scipy.ndimage.morphology import generate_binary_structure, iterate_structure
+import caiman.utils.utils
from ...paths import caiman_datadir
from .utilities import dict_compare, get_file_size
@@ -17,7 +17,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
memory_fact=1, n_processes=1, nb_patch=1, p_ssub=2, p_tsub=2,
remove_very_bad_comps=False, rf=None, stride=None,
check_nan=True, n_pixels_per_process=None,
- k=30, alpha_snmf=10e2, center_psf=False, gSig=[5, 5], gSiz=None,
+ k=30, alpha_snmf=100, center_psf=False, gSig=[5, 5], gSiz=None,
init_iter=2, method_init='greedy_roi', min_corr=.85,
min_pnr=20, gnb=1, normalize_init=True, options_local_NMF=None,
ring_size_factor=1.5, rolling_length=100, rolling_sum=True,
@@ -35,7 +35,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
sniper_mode=False, test_both=False, thresh_CNN_noisy=0.5,
thresh_fitness_delta=-50, thresh_fitness_raw=None, thresh_overlap=0.5,
update_freq=200, update_num_comps=True, use_dense=True, use_peak_max=True,
- only_init_patch=True, var_name_hdf5='mov', params_dict={},
+ only_init_patch=True, var_name_hdf5='mov', max_merge_area=None, params_dict={},
):
"""Class for setting the processing parameters. All parameters for CNMF, online-CNMF, quality testing,
and motion correction can be set here and then used in the various processing pipeline steps.
@@ -67,6 +67,18 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
var_name_hdf5: str, default: 'mov'
if loading from hdf5 name of the variable to load
+ caiman_version: str
+ version of CaImAn being used
+
+ last_commit: str
+ hash of last commit in the caiman repo
+
+ mmap_F: list[str]
+ paths to F-order memory mapped files after motion correction
+
+ mmap_C: str
+ path to C-order memory mapped file after motion correction
+
PATCH PARAMS (CNMFParams.patch)######
rf: int or None, default: None
@@ -154,6 +166,24 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
K: int, default: 30
number of components to be found (per patch or whole FOV depending on whether rf=None)
+ SC_kernel: {'heat', 'cos', binary'}, default: 'heat'
+ kernel for graph affinity matrix
+
+ SC_sigma: float, default: 1
+ variance for SC kernel
+
+ SC_thr: float, default: 0,
+ threshold for affinity matrix
+
+ SC_normalize: bool, default: True
+ standardize entries prior to computing the affinity matrix
+
+ SC_use_NN: bool, default: False
+ sparsify affinity matrix by using only nearest neighbors
+
+ SC_nnn: int, default: 20
+ number of nearest neighbors to use
+
gSig: [int, int], default: [5, 5]
radius of average neurons (in pixels)
@@ -172,6 +202,9 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
nb: int, default: 1
number of background components
+ lambda_gnmf: float, default: 1.
+ regularization weight for graph NMF
+
maxIter: int, default: 5
number of HALS iterations during initialization
@@ -208,7 +241,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
max_iter_snmf : int, default: 500
maximum number of iterations for sparse NMF initialization
- alpha_snmf: float, default: 10e2
+ alpha_snmf: float, default: 100
sparse NMF sparsity regularization weight
sigma_smooth_snmf : (float, float, float), default: (.5,.5,.5)
@@ -307,6 +340,9 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
lags: int, default: 5
number of autocovariance lags to be considered for time constant estimation
+ optimize_g: bool, default: False
+ flag for optimizing time constants
+
fudge_factor: float (close but smaller than 1) default: .96
bias correction factor for discrete time constants
@@ -332,6 +368,12 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
thr: float, default: 0.8
Trace correlation threshold for merging two components.
+ merge_parallel: bool, default: False
+ Perform merging in parallel
+
+ max_merge_area: int or None, default: None
+ maximum area (in pixels) of merged components, used to determine whether to merge components during fitting process
+
QUALITY EVALUATION PARAMETERS (CNMFParams.quality)###########
min_SNR: float, default: 2.5
@@ -532,7 +574,11 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
'fr': fr,
'decay_time': decay_time,
'dxy': dxy,
- 'var_name_hdf5': var_name_hdf5
+ 'var_name_hdf5': var_name_hdf5,
+ 'caiman_version': '1.5.2',
+ 'last_commit': None,
+ 'mmap_F': None,
+ 'mmap_C': None
}
self.patch = {
@@ -569,6 +615,14 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
self.init = {
'K': k, # number of components,
+ 'SC_kernel': 'heat', # kernel for graph affinity matrix
+ 'SC_sigma' : 1, # std for SC kernel
+ 'SC_thr': 0, # threshold for affinity matrix
+ 'SC_normalize': True, # standardize entries prior to
+ # computing affinity matrix
+ 'SC_use_NN': False, # sparsify affinity matrix by using
+ # only nearest neighbors
+ 'SC_nnn': 20, # number of nearest neighbors to use
'alpha_snmf': alpha_snmf,
'center_psf': center_psf,
'gSig': gSig,
@@ -576,6 +630,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
'gSiz': gSiz,
'init_iter': init_iter,
'kernel': None, # user specified template for greedyROI
+ 'lambda_gnmf' :1, # regularization weight for graph NMF
'maxIter': 5, # number of HALS iterations
'max_iter_snmf': 500,
'method_init': method_init, # can be greedy_roi, greedy_pnr sparse_nmf, local_NMF
@@ -616,8 +671,8 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
'normalize_yyt_one': True,
'nrgthr': 0.9999, # Energy threshold
'num_blocks_per_run_spat': num_blocks_per_run_spat, # number of process to parallelize residual computation ** DECREASE IF MEMORY ISSUES
- 'se': None, # Morphological closing structuring element
- 'ss': None, # Binary element for determining connectivity
+ 'se': np.ones((3, 3), dtype='uint8'), # Morphological closing structuring element
+ 'ss': np.ones((3, 3), dtype='uint8'), # Binary element for determining connectivity
'thr_method': 'nrg', # Method of thresholding ('max' or 'nrg')
# whether to update the background components in the spatial phase
'update_background_components': update_background_components,
@@ -634,6 +689,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
'fudge_factor': .96,
# number of autocovariance lags to be considered for time constant estimation
'lags': 5,
+ 'optimize_g': False, # flag for optimizing time constants
'memory_efficient': False,
# method for solving the constrained deconvolution problem ('oasis','cvx' or 'cvxpy')
# if method cvxpy, primary and secondary (if problem unfeasible for approx
@@ -652,6 +708,8 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
self.merging = {
'do_merge': do_merge,
'merge_thr': merge_thresh,
+ 'merge_parallel': False,
+ 'max_merge_area': max_merge_area
}
self.quality = {
@@ -729,6 +787,7 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
}
self.change_params(params_dict)
+ self.data['last_commit'] = '-'.join(caiman.utils.utils.get_caiman_version())
if self.data['dims'] is None and self.data['fnames'] is not None:
self.data['dims'] = get_file_size(self.data['fnames'], var_name_hdf5=self.data['var_name_hdf5'])[0]
if self.data['fnames'] is not None:
@@ -753,7 +812,8 @@ def __init__(self, fnames=None, dims=None, dxy=(1, 1),
if self.init['gSig'] is None:
self.init['gSig'] = [-1, -1]
if self.init['gSiz'] is None:
- self.init['gSiz'] = [2*gs + 1 for gs in self.init['gSig']]
+ self.init['gSiz'] = [2*gs + 3 for gs in self.init['gSig']]
+ self.init['gSiz'] = [gz if gz % 2 else gz + 1 for gz in self.init['gSiz']]
if gnb <= 0:
logging.warning("gnb={0}, hence setting keys nb_patch and low_rank_background ".format(gnb) +
diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py
index 12a175b14..c06e62daa 100644
--- a/caiman/source_extraction/cnmf/spatial.py
+++ b/caiman/source_extraction/cnmf/spatial.py
@@ -31,6 +31,7 @@
from sklearn.decomposition import NMF
import tempfile
import time
+import psutil
from typing import List
from ...mmapping import load_memmap, parallel_dot_product
@@ -43,8 +44,6 @@ def basis_denoising(y, c, boh, sn, id2_, px):
else:
return (None, None, None)
return a, px, id2_
-#%% update_spatial_components (in parallel)
-
def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None,
min_size=3, max_size=8, dist=3,
@@ -193,6 +192,14 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None,
ff = ff[0]
A_in = csc_column_remove(A_in, list(ff))
C = np.delete(C, list(ff), 0)
+
+ # remove empty components
+ ff = np.where(np.sum(C, axis=1)==0)
+ if np.size(ff) > 0:
+ logging.info("Eliminating empty components: {}".format(ff))
+ ff = ff[0]
+ A_in = csc_column_remove(A_in, list(ff))
+ C = np.delete(C, list(ff), 0)
if normalize_yyt_one and C is not None:
C = np.array(C)
@@ -220,7 +227,6 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None,
if i + n_pixels_per_process < np.prod(dims):
pixel_groups.append([Y_name, C_name, sn, ind2_[(i + n_pixels_per_process):np.prod(dims)], list(
range(i + n_pixels_per_process, np.prod(dims))), method_ls, cct])
- #A_ = np.zeros((d, nr + np.size(f, 0))) # init A_
#A_ = scipy.sparse.lil_matrix((d, nr + np.size(f, 0)))
if dview is not None:
if 'multiprocessing' in str(type(dview)):
@@ -232,9 +238,9 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None,
dview.results.clear()
else:
parallel_result = list(map(regression_ipyparallel, pixel_groups))
- data = []
- rows = []
- cols = []
+ data:List = []
+ rows:List = []
+ cols:List = []
for chunk in parallel_result:
for pars in chunk:
px, idxs_, a = pars
@@ -266,18 +272,18 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None,
f = np.delete(f, background_ff, 0)
b_in = np.delete(b_in, background_ff, 1)
- A_ = A_[:, :nr]
- A_ = coo_matrix(A_)
- logging.info("Computing residuals")
-
- if 'memmap' in str(type(Y)):
- Y_resf = parallel_dot_product(Y, f.T, dview=dview, block_size=block_size_spat, num_blocks_per_run=num_blocks_per_run_spat) - \
- A_.dot(C[:nr].dot(f.T))
- else:
- # Y*f' - A*(C*f')
- Y_resf = np.dot(Y, f.T) - A_.dot(C[:nr].dot(f.T))
-
+ A_ = A_[:, :nr]
if update_background_components:
+ A_ = csr_matrix(A_)
+ logging.info("Computing residuals")
+ if 'memmap' in str(type(Y)):
+ bl_siz1 = Y.shape[0] // (num_blocks_per_run_spat - 1)
+ bl_siz2 = psutil.virtual_memory().available // (4*Y.shape[-1]*(num_blocks_per_run_spat + 1))
+ Y_resf = parallel_dot_product(Y, f.T, dview=dview, block_size=min(bl_siz1, bl_siz2), num_blocks_per_run=num_blocks_per_run_spat) - \
+ A_.dot(C[:nr].dot(f.T))
+ else:
+ # Y*f' - A*(C*f')
+ Y_resf = np.dot(Y, f.T) - A_.dot(C[:nr].dot(f.T))
if b_in is None:
# update baseline based on residual
@@ -305,10 +311,8 @@ def update_spatial_components(Y, C=None, f=None, A_in=None, sn=None, dims=None,
except:
raise Exception("Failed to delete: " + folder)
- return A_, b, C, f
+ return csc_matrix(A_), b, C, f
-
-#%%
def HALS4shape_bckgrnd(Y_resf, B, F, ind_B, iters=5):
K = B.shape[-1]
U = Y_resf.T
@@ -323,7 +327,6 @@ def HALS4shape_bckgrnd(Y_resf, B, F, ind_B, iters=5):
return B
-# %%lars_regression_noise_ipyparallel
def regression_ipyparallel(pars):
"""update spatial footprints and background through Basis Pursuit Denoising
@@ -386,7 +389,6 @@ def regression_ipyparallel(pars):
_, T = np.shape(C) # 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], :]
idx_only_neurons = idxs_C[idx_px_from_0]
@@ -407,7 +409,11 @@ def regression_ipyparallel(pars):
elif method_least_square == 'lasso_lars': # lasso lars function from scikit learn
lambda_lasso = 0 if np.size(cct_) == 0 else \
.5 * noise_sn[px] * np.sqrt(np.max(cct_)) / T
- clf = linear_model.LassoLars(alpha=lambda_lasso, positive=True, fit_intercept=True)
+ clf = linear_model.LassoLars(alpha=lambda_lasso, positive=True,
+ fit_intercept=True)
+# clf = linear_model.Lasso(alpha=lambda_lasso, positive=True,
+# fit_intercept=True, normalize=True,
+# selection='random')
a_lrs = clf.fit(np.array(c.T), np.ravel(y))
a = a_lrs.coef_
@@ -429,7 +435,6 @@ def regression_ipyparallel(pars):
return As
-# %%
def construct_ellipse_parallel(pars):
"""update spatial footprints and background through Basis Pursuit Denoising
@@ -452,10 +457,6 @@ def construct_ellipse_parallel(pars):
# search indexes for each component
return np.sqrt(np.sum([old_div((dist_cm * V[:, k]) ** 2, dkk[k]) for k in range(len(dkk))], 0)) <= dist
-
-# %% threshold_components
-
-
def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgthr=0.9999, extract_cc=True,
se=None, ss=None, dview=None):
"""
@@ -527,9 +528,9 @@ def threshold_components(A, dims, medw=None, thr_method='max', maxthr=0.1, nrgth
res = list(map(threshold_components_parallel, pars))
res.sort(key=lambda x: x[1])
- indices = []
+ indices:List = []
indptr = [0]
- data = []
+ data:List = []
for r in res:
At, i = r
indptr.append(indptr[-1]+At.indptr[-1])
@@ -632,10 +633,6 @@ def threshold_components_parallel(pars):
return csr_matrix(Ath2), i
-
-# %%
-
-
def nnls_L0(X, Yp, noise):
"""
Nonnegative least square with L0 penalty
@@ -676,7 +673,9 @@ def nnls_L0(X, Yp, noise):
else:
W_lam[eliminate[np.argmin(np.array(eliminate)[:, 1])][0]] = 0
-# %% auxiliary functions
+#####
+# auxiliary functions
+
def calcAvec(new, dQ, W, lambda_, active_set, M, positive):
"""
calculate the vector to travel along
@@ -917,7 +916,7 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8,
elif method == 'dilate':
indptr = [0]
- indices = []
+ indices:List = []
data = []
if dview is None:
for i in range(nr):
@@ -964,10 +963,8 @@ def determine_search_location(A, dims, method='ellipse', min_size=3, max_size=8,
raise Exception('Not implemented')
dist_indicator = True * np.ones((d, nr))
-
-
return csc_matrix(dist_indicator)
-#%%
+
def construct_dilate_parallel(pars):
"""
"""
@@ -988,7 +985,7 @@ def construct_dilate_parallel(pars):
# search indexes for each component
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)
@@ -1100,9 +1097,6 @@ def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size,
return ind2_, nr, C, f, b, A_in
-
-
-#%%
def creatememmap(Y, Cf, dview):
"""memmap the C and Y objects in parallel
@@ -1127,7 +1121,7 @@ def creatememmap(Y, Cf, dview):
"""
if os.environ.get('SLURM_SUBMIT_DIR') is not None:
tmpf = os.environ.get('SLURM_SUBMIT_DIR')
- print(('cluster temporary folder:' + tmpf))
+ print(f'cluster temporary folder: {tmpf}')
folder = tempfile.mkdtemp(dir=tmpf)
else:
folder = tempfile.mkdtemp()
diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py
index 29d0c4f4a..dafc83306 100644
--- a/caiman/source_extraction/cnmf/temporal.py
+++ b/caiman/source_extraction/cnmf/temporal.py
@@ -9,10 +9,11 @@
from builtins import map
from builtins import range
import logging
-from scipy.sparse import spdiags, diags, coo_matrix # ,csgraph
+from scipy.sparse import spdiags, diags, coo_matrix, csc_matrix # ,csgraph
import scipy
import numpy as np
import platform
+import psutil
from .deconvolution import constrained_foopsi
from .utilities import update_order_greedy
import sys
@@ -194,24 +195,25 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N
A = scipy.sparse.hstack((A, b)).tocsc()
S = np.zeros(np.shape(Cin))
Cin = np.vstack((Cin, fin))
- C = Cin.copy()
+ #C = Cin.copy()
nA = np.ravel(A.power(2).sum(axis=0))
logging.info('Generating residuals')
# dview_res = None if block_size >= 500 else dview
if 'memmap' in str(type(Y)):
- YA = parallel_dot_product(Y, A, dview=dview, block_size=block_size_temp,
- transpose=True, num_blocks_per_run=num_blocks_per_run_temp) * diags(1. / nA)
+ bl_siz1 = d // (np.maximum(num_blocks_per_run_temp - 1, 1))
+ bl_siz2 = int(psutil.virtual_memory().available/(num_blocks_per_run_temp + 1) - 4*A.nnz) // int(4*T)
+ # block_size_temp
+ YA = parallel_dot_product(Y, A.tocsr(), dview=dview, block_size=min(bl_siz1, bl_siz2),
+ transpose=True, num_blocks_per_run=num_blocks_per_run_temp) * diags(1. / nA);
else:
- # YA = (A.T.dot(Y).T) * spdiags(old_div(1., nA), 0, nr + nb, nr + nb)
YA = (A.T.dot(Y).T) * diags(1. / nA)
- # AA = ((A.T.dot(A)) * spdiags(old_div(1., nA), 0, nr + nb, nr + nb)).tocsr()
AA = ((A.T.dot(A)) * diags(1. / nA)).tocsr()
YrA = YA - AA.T.dot(Cin).T
# creating the patch of components to be computed in parrallel
parrllcomp, len_parrllcomp = update_order_greedy(AA[:nr, :][:, :nr])
logging.info("entering the deconvolution ")
- C, S, bl, YrA, c1, sn, g, lam = update_iteration(parrllcomp, len_parrllcomp, nb, C, S, bl, nr,
+ C, S, bl, YrA, c1, sn, g, lam = update_iteration(parrllcomp, len_parrllcomp, nb, Cin, 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
@@ -235,7 +237,7 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N
nr = nr - (len(ff) - len(background_ff))
b = A[:, nr:].toarray()
- A = coo_matrix(A[:, :nr])
+ A = csc_matrix(A[:, :nr])
f = C[nr:, :]
C = C[:nr, :]
YrA = np.array(YrA[:, :nr]).T
diff --git a/caiman/source_extraction/cnmf/utilities.py b/caiman/source_extraction/cnmf/utilities.py
index 155f1be5e..ec3735cdb 100644
--- a/caiman/source_extraction/cnmf/utilities.py
+++ b/caiman/source_extraction/cnmf/utilities.py
@@ -38,8 +38,8 @@
from .initialization import greedyROI
from ...base.rois import com
-
from ...mmapping import parallel_dot_product, load_memmap
+from ...cluster import extract_patch_coordinates
from ...utils.stats import df_percentile
@@ -652,7 +652,6 @@ def manually_refine_components(Y, xxx_todo_changeme, A, C, Cn, thr=0.9, display_
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)
-
y3_res = np.reshape(y2_res, (dy_sz, dx_sz, T), order='F')
a__, c__, center__, b_in__, f_in__ = greedyROI(
y3_res, nr=1, gSig=[np.floor(old_div(dx_sz, 2)), np.floor(old_div(dy_sz, 2))], gSiz=[dx_sz, dy_sz])
@@ -986,13 +985,13 @@ def get_file_size(file_name, var_name_hdf5='mov'):
dims = [0, 0]
try:
T = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))
- dims[0] = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
- dims[1] = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))
+ dims[1] = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
+ dims[0] = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))
except():
print('Roll back to opencv 2')
T = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_COUNT))
- dims[0] = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_WIDTH))
- dims[1] = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT))
+ dims[1] = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_WIDTH))
+ dims[0] = int(cap.get(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT))
elif extension == '.mmap':
filename = os.path.split(file_name)[-1]
Yr, dims, T = load_memmap(os.path.join(
@@ -1048,3 +1047,146 @@ def get_file_size(file_name, var_name_hdf5='mov'):
else:
raise Exception('Unknown input type')
return dims, T
+
+
+def fast_graph_Laplacian(mmap_file, dims, max_radius=10, kernel='heat',
+ dview=None, sigma=1, thr=0.05, p=10, normalize=True,
+ use_NN=False, rf=None, strides=None):
+ """ Computes an approximate affinity maps and its graph Laplacian for all
+ pixels. For each pixel it restricts its attention to a given radius around
+ it.
+ Args:
+ mmap_file: str
+ Memory mapped file in pixel first order
+
+ max_radius: float
+ Maximum radius around each pixel
+
+ kernel: str {'heat', 'binary', 'cos'}
+ type of kernel
+
+ dview: dview object
+ multiprocessing or ipyparallel object for parallelization
+
+ sigma: float
+ standard deviation of Gaussian (heat) kernel
+
+ thr: float
+ threshold for affinity matrix
+
+ p: int
+ number of neighbors
+
+ normalize: bool
+ normalize vectors before computing affinity
+
+ use_NN: bool
+ use only p nearest neighbors
+
+ Returns:
+ W: scipy.sparse.csr_matrix
+ Graph affinity matrix
+
+ D: scipy.sparse.spdiags
+ Diagonal of affinity matrix
+
+ L: scipy.sparse.csr_matrix
+ Graph Laplacian matrix
+ """
+ Np = np.prod(np.array(dims))
+ if rf is None:
+ pars = []
+ for i in range(Np):
+ pars.append([i, mmap_file, dims, max_radius, kernel, sigma, thr,
+ p, normalize, use_NN])
+ if dview is None:
+ res = list(map(fast_graph_Laplacian_pixel, pars))
+ else:
+ res = dview.map(fast_graph_Laplacian_pixel, pars, chunksize=128)
+ indptr = np.cumsum(np.array([0] + [len(r[0]) for r in res]))
+ indeces = [item for sublist in res for item in sublist[0]]
+ data = [item for sublist in res for item in sublist[1]]
+ W = scipy.sparse.csr_matrix((data, indeces, indptr), shape=[Np, Np])
+ D = scipy.sparse.spdiags(W.sum(0), 0, Np, Np)
+ L = D - W
+ else:
+ indices, _ = extract_patch_coordinates(dims, rf, strides)
+ pars = []
+ for i in range(len(indices)):
+ pars.append([mmap_file, indices[i], kernel, sigma, thr, p,
+ normalize, use_NN])
+ if dview is None:
+ res = list(map(fast_graph_Laplacian_patches, pars))
+ else:
+ res = dview.map(fast_graph_Laplacian_patches, pars)
+ W = res
+ D = [scipy.sparse.spdiags(w.sum(0), 0, w.shape[0], w.shape[0]) for w in W]
+ L = [d - w for (d, w) in zip(W, D)]
+ return W, D, L
+
+
+def fast_graph_Laplacian_patches(pars):
+ """ Computes the full graph affinity matrix on a patch. See
+ fast_graph_Laplacian above for definition of arguments.
+ """
+ mmap_file, indices, kernel, sigma, thr, p, normalize, use_NN = pars
+ if type(mmap_file) not in {'str', 'list'}:
+ Yind = mmap_file
+ else:
+ Y = load_memmap(mmap_file)[0]
+ Yind = np.array(Y[indices])
+ if normalize:
+ Yind -= Yind.mean(1)[:, np.newaxis]
+ Yind /= np.sqrt((Yind**2).sum(1)[:, np.newaxis])
+ yf = np.ones((Yind.shape[0], 1))
+ else:
+ yf = (Yind**2).sum(1)[:, np.newaxis]
+ yyt = Yind.dot(Yind.T)
+ W = np.exp(-(yf + yf.T - 2*yyt)/sigma) if kernel.lower() == 'heat' else yyt
+ W[W\\n\"+\n", + " \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n", + " \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n", + " \"
\\n\"+\n", + " \"\\n\"+\n",
+ " \"from bokeh.resources import INLINE\\n\"+\n",
+ " \"output_notebook(resources=INLINE)\\n\"+\n",
+ " \"\\n\"+\n",
+ " \"\\n\"+\n \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n \"
\\n\"+\n \"\\n\"+\n \"from bokeh.resources import INLINE\\n\"+\n \"output_notebook(resources=INLINE)\\n\"+\n \"\\n\"+\n \"