In [ ]:
#%%
try:
get_ipython().magic(u'load_ext autoreload')
get_ipython().magic(u'autoreload 2')
print(1)
except:
print('NOT IPYTHON')
from ipyparallel import Client
import logging
import matplotlib.pyplot as plt
import numpy as np
import os
import psutil
from scipy.ndimage.filters import gaussian_filter
import sys
import caiman as cm
from caiman.utils.visualization import nb_view_patches3d
import caiman.source_extraction.cnmf as cnmf
from caiman.components_evaluation import evaluate_components, estimate_components_quality_auto
from caiman.cluster import setup_cluster
from caiman.paths import caiman_datadir
import bokeh.plotting as bpl
bpl.output_notebook()
logging.basicConfig(format=
"%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s",
# filename="/tmp/caiman.log",
level=logging.DEBUG)
In [ ]:
# stop the cluster if one exists
n_processes = psutil.cpu_count()
print('using ' + str(n_processes) + ' processes')
print("Stopping cluster to avoid unnencessary use of memory....")
sys.stdout.flush()
cm.stop_server()
Define a function to create some toy data
In [ ]:
def gen_data(p=1, noise=1., T=256, framerate=30, firerate=2., plot=False):
if p == 2:
gamma = np.array([1.5, -.55])
elif p == 1:
gamma = np.array([.9])
else:
raise
dims = (30, 40, 50) # size of image
sig = (2, 2, 2) # neurons size
bkgrd = 10
N = 20 # number of neurons
np.random.seed(7)
centers = np.asarray([[np.random.randint(5, x - 5)
for x in dims] for i in range(N)])
Yr = np.zeros(dims + (T,), dtype=np.float32)
trueSpikes = np.random.rand(N, T) < firerate / float(framerate)
trueSpikes[:, 0] = 0
truth = trueSpikes.astype(np.float32)
for i in range(2, T):
if p == 2:
truth[:, i] += gamma[0] * truth[:, i - 1] + gamma[1] * truth[:, i - 2]
else:
truth[:, i] += gamma[0] * truth[:, i - 1]
for i in range(N):
Yr[centers[i, 0], centers[i, 1], centers[i, 2]] = truth[i]
tmp = np.zeros(dims)
tmp[15, 20, 25] = 1.
z = np.linalg.norm(gaussian_filter(tmp, sig).ravel())
Yr = bkgrd + noise * np.random.randn(*(dims + (T,))) + 10 * gaussian_filter(Yr, sig + (0,)) / z
d1, d2, d3, T = Yr.shape
Yr = np.reshape(Yr, (d1 * d2 * d3, T), order='F').astype(np.float32)
if plot:
Y = np.reshape(Yr, (d1, d2, d3, T), order='F')
plt.figure(figsize=(15, 3))
plt.plot(truth.T)
plt.figure(figsize=(15, 3))
for c in centers:
plt.plot(Y[c[0], c[1], c[2]])
plt.figure(figsize=(15, 4))
plt.subplot(131)
plt.scatter(*centers.T[::-1], c='g')
plt.imshow(Y.max(0).max(-1), cmap='hot')
plt.title('Max.proj. x & t')
plt.subplot(132)
plt.scatter(*centers.T[[2, 0, 1]], c='g')
plt.imshow(Y.max(1).max(-1), cmap='hot')
plt.title('Max.proj. y & t')
plt.subplot(133)
plt.scatter(*centers.T[[1, 0, 2]], c='g')
plt.imshow(Y.max(2).max(-1), cmap='hot')
plt.title('Max.proj. z & t')
plt.show()
return Yr, truth, trueSpikes, centers, dims
In [ ]:
plt.close('all')
#%% SAVING TIFF FILE ON A SINGLE MEMORY MAPPABLE FILE
demo_filename = os.path.join(caiman_datadir(), 'example_movies', 'demoMovie3D.tif')
try:
fname_new = cm.save_memmap([demo_filename], base_name='Yr', is_3D=True, order='C')
except: # %% create 3d tiff file if not yet existent
from skimage.external.tifffile import imsave
Yr, truth, trueSpikes, centers, dims = gen_data(p=2)
data = np.transpose(Yr.reshape(dims + (-1,), order='F'), [3, 0, 1, 2])
imsave(demo_filename, data)
fname_new = cm.save_memmap([demo_filename], base_name='Yr', is_3D=True, order='C')
print(fname_new)
Load memory mapped file and show a max-projection of the correlation image
In [ ]:
Yr, dims, T = cm.load_memmap(fname_new)
Y = np.reshape(Yr, dims + (T,), order='F')
Cn = cm.local_correlations(Y)
plt.imshow(Cn.max(0) if len(Cn.shape) == 3 else Cn, cmap='gray',
vmin=np.percentile(Cn, 1), vmax=np.percentile(Cn, 99))
plt.show()
In [ ]:
# set parameters
K = 20 # number of neurons expected per patch
gSig = [2, 2, 2] # expected half size of neurons
merge_thresh = 0.8 # merging threshold, max correlation allowed
p = 2 # order of the autoregressive system
In [ ]:
%%capture
# START CLUSTER
c, dview, n_processes = setup_cluster(
backend='local', n_processes=None, single_thread=False)
Initialize CNMF object
In [ ]:
# INIT
cnm = cnmf.CNMF(n_processes, method_init='greedy_roi', k=K, gSig=gSig, merge_thresh=merge_thresh,
p=p, dview=dview, Ain=None, method_deconvolution='oasis')
In [ ]:
%%capture
# FIT
images = np.reshape(Yr.T, [T] + list(dims), order='F') # reshape data in Python format (T x X x Y x Z)
cnm = cnm.fit(images)
In [ ]:
cnm.estimates.nb_view_components_3d(image_type='mean', dims=dims)
In [ ]:
%%capture
rf = (15, 15, 15) # half-size of the patches in pixels. rf=25, patches are 50x50
stride = (10, 10, 10) # amounpl.it of overlap between the patches in pixels
K = 12 # number of neurons expected per patch
gSig = [2, 2, 2] # expected half size of neurons
merge_thresh = 0.8 # merging threshold, max correlation allowed
p = 2 # order of the autoregressive system
save_results = False
#%% RUN ALGORITHM ON PATCHES
init_method = 'greedy_roi'
alpha_snmf = None # 10e2 # this controls sparsity
cnm = cnmf.CNMF(n_processes, k=K, gSig=gSig, merge_thresh=0.8, p=p, dview=dview, Ain=None, rf=rf, stride=stride, memory_fact=1,
method_init=init_method, alpha_snmf=alpha_snmf, only_init_patch=True, gnb=1, method_deconvolution='oasis')
cnm = cnm.fit(images)
print(('Number of components:' + str(cnm.estimates.A.shape[-1])))
In [ ]:
#%% COMPONENT EVALUATION
# the components are evaluated in two ways:
# a) the shape of each component must be correlated with the data
# b) a minimum peak SNR is required over the length of a transient
fr = 10 # approx final rate (after eventual downsampling )
decay_time = 1. # length of typical transient in seconds
use_cnn = False # CNN classifier is designed for 2d (real) data
min_SNR = 3 # accept components with that peak-SNR or higher
rval_thr = 0.7 # accept components iwth speace correlation threshold or higher
cnm.params.change_params(params_dict={'fr': fr,
'decay_time': decay_time,
'min_SNR': min_SNR,
'rval_thr': rval_thr,
'use_cnn': use_cnn})
cnm.estimates.evaluate_components(images, cnm.params, dview=dview)
In [ ]:
print(('Keeping ' + str(len(cnm.estimates.idx_components)) +
' and discarding ' + str(len(cnm.estimates.idx_components_bad))))
In [ ]:
%%capture
#%% RE-RUN seeded CNMF on accepted components
cnm.params.set('temporal', {'p': p})
cnm2 = cnm.refit(images)
In [ ]:
# view components per layer
cnm2.estimates.nb_view_components_3d(image_type='corr', dims=dims, Yr=Yr)
In [ ]:
# STOP CLUSTER
cm.stop_server(dview=dview)