The demo shows how to perform online analysis on one photon data using a Ring CNN for extracting the background followed by processing using the OnACID algorithm. The algorithm relies on the usage a GPU to efficiently estimate and apply the background model so it is recommended to have access to a GPU when running this notebook.
In [ ]:
get_ipython().magic('load_ext autoreload')
get_ipython().magic('autoreload 2')
from IPython.display import display, clear_output
import glob
import logging
import numpy as np
import os
import cv2
logging.basicConfig(format=
"%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s",
# filename="/tmp/caiman.log",
level=logging.WARNING)
import caiman as cm
from caiman.source_extraction import cnmf as cnmf
from caiman.paths import caiman_datadir
from caiman.utils.utils import download_demo
import matplotlib.pyplot as plt
import bokeh.plotting as bpl
bpl.output_notebook()
In [ ]:
fnames=download_demo('blood_vessel_10Hz.mat')
In [ ]:
reuse_model = False # set to True to re-use an existing ring model
path_to_model = None # specify a pre-trained model here if needed
gSig = (7, 7) # expected half size of neurons
gnb = 2 # number of background components for OnACID
init_batch = 500 # number of frames for initialization and training
params_dict = {'fnames': fnames,
'var_name_hdf5': 'Y', # name of variable inside mat file where the data is stored
'fr': 10, # frame rate (Hz)
'decay_time': 0.5, # approximate length of transient event in seconds
'gSig': gSig,
'p': 0, # order of AR indicator dynamics
'ring_CNN': True, # SET TO TRUE TO USE RING CNN
'min_SNR': 2.65, # minimum SNR for accepting new components
'SNR_lowest': 0.75, # reject components with SNR below this value
'use_cnn': False, # do not use CNN based test for components
'use_ecc': True, # test eccentricity
'max_ecc': 2.625, # reject components with eccentricity above this value
'rval_thr': 0.70, # correlation threshold for new component inclusion
'rval_lowest': 0.25, # reject components with corr below that value
'ds_factor': 1, # spatial downsampling factor (increases speed but may lose some fine structure)
'nb': gnb,
'motion_correct': False, # Flag for motion correction
'init_batch': init_batch, # number of frames for initialization (presumably from the first file)
'init_method': 'bare',
'normalize': False,
'expected_comps': 1100, # maximum number of expected components used for memory pre-allocation (exaggerate here)
'sniper_mode': False, # flag using a CNN to detect new neurons (o/w space correlation is used)
'dist_shape_update' : True, # flag for updating shapes in a distributed way
'min_num_trial': 5, # number of candidate components per frame
'epochs': 3, # number of total passes over the data
'stop_detection': True, # Run a last epoch without detecting new neurons
'K': 50, # initial number of components
'lr': 6e-4,
'lr_scheduler': [0.9, 6000, 10000],
'pct': 0.01,
'path_to_model': path_to_model, # where the ring CNN model is saved/loaded
'reuse_model': reuse_model # flag for re-using a ring CNN model
}
opts = cnmf.params.CNMFParams(params_dict=params_dict)
The first initbatch frames are used for training the ring-CNN model. Once the model is trained the background is subtracted and the different is used for initialization purposes. The initialization method chosen here bare will only search for a small number of neurons and is mostly used to initialize the background components. Initialization with the full CNMF can also be used by choosing cnmf.
We first create an OnACID object located in the module online_cnmf and we pass the parameters similarly to the case of batch processing. We then run the algorithm using the fit_online method. We then save the results inside
the folder where the Ring_CNN model is saved.
In [ ]:
run_onacid = True
if run_onacid:
cnm = cnmf.online_cnmf.OnACID(params=opts)
cnm.fit_online()
fld_name = os.path.dirname(cnm.params.ring_CNN['path_to_model'])
res_name_nm = os.path.join(fld_name, 'onacid_results_nm.hdf5')
cnm.save(res_name_nm) # save initial results (without any postprocessing)
else:
fld_name = os.path.dirname(path_to_model)
res_name = os.path.join(fld_name, 'onacid_results.hdf5')
cnm = cnmf.online_cnmf.load_OnlineCNMF(res_name)
cnm.params.data['fnames'] = fnames
In [ ]:
ds = 10 # plot every ds frames to make more manageable figures
init_batch = 500
dims, T = cnmf.utilities.get_file_size(fnames, var_name_hdf5='Y')
T = np.array(T).sum()
n_epochs = cnm.params.online['epochs']
T_detect = 1e3*np.hstack((np.zeros(init_batch), cnm.t_detect))
T_shapes = 1e3*np.hstack((np.zeros(init_batch), cnm.t_shapes))
T_online = 1e3*np.hstack((np.zeros(init_batch), cnm.t_online)) - T_detect - T_shapes
plt.figure()
plt.stackplot(np.arange(len(T_detect))[::ds], T_online[::ds], T_detect[::ds], T_shapes[::ds],
colors=['tab:red', 'tab:purple', 'tab:brown'])
plt.legend(labels=['process', 'detect', 'shapes'], loc=2)
plt.title('Processing time allocation')
plt.xlabel('Frame #')
plt.ylabel('Processing time [ms]')
max_val = 80
plt.ylim([0, max_val]);
plt.plot([init_batch, init_batch], [0, max_val], '--k')
for i in range(n_epochs - 1):
plt.plot([(i+1)*T, (i+1)*T], [0, max_val], '--k')
plt.xlim([0, n_epochs*T]);
plt.savefig(os.path.join(fld_name, 'time_per_frame_ds.pdf'), bbox_inches='tight', pad_inches=0)
In [ ]:
init_batch = 500
plt.figure()
tc_init = cnm.t_init*np.ones(T*n_epochs)
ds = 10
#tc_mot = np.hstack((np.zeros(init_batch), np.cumsum(T_motion)/1000))
tc_prc = np.cumsum(T_online)/1000#np.hstack((np.zeros(init_batch), ))
tc_det = np.cumsum(T_detect)/1000#np.hstack((np.zeros(init_batch), ))
tc_shp = np.cumsum(T_shapes)/1000#np.hstack((np.zeros(init_batch), ))
plt.stackplot(np.arange(len(tc_init))[::ds], tc_init[::ds], tc_prc[::ds], tc_det[::ds], tc_shp[::ds],
colors=['g', 'tab:red', 'tab:purple', 'tab:brown'])
plt.legend(labels=['initialize', 'process', 'detect', 'shapes'], loc=2)
plt.title('Processing time allocation')
plt.xlabel('Frame #')
plt.ylabel('Processing time [s]')
max_val = (tc_prc[-1] + tc_det[-1] + tc_shp[-1] + cnm.t_init)*1.05
for i in range(n_epochs - 1):
plt.plot([(i+1)*T, (i+1)*T], [0, max_val], '--k')
plt.xlim([0, n_epochs*T]);
plt.ylim([0, max_val])
plt.savefig(os.path.join(fld_name, 'time_cumulative_ds.pdf'), bbox_inches='tight', pad_inches=0)
In [ ]:
print('Cost of estimating model and running first epoch: {:.2f}s'.format(tc_prc[T] + tc_det[T] + tc_shp[T] + tc_init[T]))
In [ ]:
# first compute background summary images
images = cm.load(fnames, var_name_hdf5='Y', subindices=slice(None, None, 2))
cn_filter, pnr = cm.summary_images.correlation_pnr(images, gSig=3, swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile
In [ ]:
plt.figure(figsize=(15, 7))
plt.subplot(1,2,1); plt.imshow(cn_filter); plt.colorbar()
plt.subplot(1,2,2); plt.imshow(pnr); plt.colorbar()
In [ ]:
cnm.estimates.plot_contours_nb(img=cn_filter, idx=cnm.estimates.idx_components, line_color='white', thr=0.3)
Now inspect the components extracted by OnACID. Note that if single pass was used then several components would be non-zero only for the part of the time interval indicating that they were detected online by OnACID.
Note that if you get data rate error you can start Jupyter notebooks using: 'jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10'
In [ ]:
cnm.estimates.nb_view_components(img=cn_filter, denoised_color='red')
In [ ]:
save_file = True
if save_file:
from caiman.utils.nn_models import create_LN_model
model_LN = create_LN_model(images, shape=opts.data['dims'] + (1,), n_channels=opts.ring_CNN['n_channels'],
width=opts.ring_CNN['width'], use_bias=opts.ring_CNN['use_bias'], gSig=gSig[0],
use_add=opts.ring_CNN['use_add'])
model_LN.load_weights(cnm.params.ring_CNN['path_to_model'])
# Load the data in batches and save them
m = []
saved_files = []
batch_length = 256
for i in range(0, T, batch_length):
images = cm.load(fnames, var_name_hdf5='Y', subindices=slice(i, i + batch_length))
images_filt = np.squeeze(model_LN.predict(np.expand_dims(images, axis=-1)))
temp_file = os.path.join(fld_name, 'pfc_back_removed_' + format(i, '05d') + '.h5')
saved_files.append(temp_file)
m = cm.movie(np.maximum(images - images_filt, 0))
m.save(temp_file)
else:
saved_files = glob.glob(os.path.join(fld_name, 'pfc_back_removed_*'))
saved_files.sort()
In [ ]:
fname_mmap = cm.save_memmap([saved_files], order='C', border_to_0=0)
Yr, dims, T = cm.load_memmap(fname_mmap)
images_mmap = Yr.T.reshape((T,) + dims, order='F')
In [ ]:
cnm.params.merging['merge_thr'] = 0.7
cnm.estimates.c1 = np.zeros(cnm.estimates.A.shape[-1])
cnm.estimates.bl = np.zeros(cnm.estimates.A.shape[-1])
cnm.estimates.neurons_sn = np.zeros(cnm.estimates.A.shape[-1])
cnm.estimates.g = None #np.ones((cnm.estimates.A.shape[-1], 1))*.9
cnm.estimates.merge_components(Yr, cnm.params)
In [ ]:
cnm.params.quality
In [ ]:
cnm.estimates.evaluate_components(imgs=images_mmap, params=cnm.params)
In [ ]:
cnm.estimates.plot_contours_nb(img=cn_filter, idx=cnm.estimates.idx_components, line_color='white')
In [ ]:
cnm.estimates.nb_view_components(idx=cnm.estimates.idx_components, img=cn_filter)
In [ ]:
cnmfe_results = download_demo('online_vs_offline.npz')
locals().update(np.load(cnmfe_results, allow_pickle=True))
A_patch_good = A_patch_good.item()
estimates_gt = cnmf.estimates.Estimates(A=A_patch_good, C=C_patch_good, dims=dims)
In [ ]:
maxthr=0.01
cnm.estimates.A_thr=None
cnm.estimates.threshold_spatial_components(maxthr=maxthr)
estimates_gt.A_thr=None
estimates_gt.threshold_spatial_components(maxthr=maxthr*10)
min_size = np.pi*(gSig[0]/1.5)**2
max_size = np.pi*(gSig[0]*1.5)**2
ntk = cnm.estimates.remove_small_large_neurons(min_size_neuro=min_size, max_size_neuro=2*max_size)
gtk = estimates_gt.remove_small_large_neurons(min_size_neuro=min_size, max_size_neuro=2*max_size)
In [ ]:
m1, m2, nm1, nm2, perf = cm.base.rois.register_ROIs(estimates_gt.A_thr[:, estimates_gt.idx_components],
cnm.estimates.A_thr[:, cnm.estimates.idx_components],
dims, align_flag=False, thresh_cost=.7, plot_results=True,
Cn=cn_filter, enclosed_thr=None)[:-1]
In [ ]:
for k, v in perf.items():
print(k + ':', '%.4f' % v, end=' ')
In [ ]:
res_name = os.path.join(fld_name, 'onacid_results.hdf5')
cnm.save(res_name)
In [ ]:
import matplotlib.lines as mlines
lp, hp = np.nanpercentile(cn_filter, [5, 98])
A_onacid = cnm.estimates.A_thr.toarray().copy()
A_onacid /= A_onacid.max(0)
A_TP = estimates_gt.A[:, m1].toarray() #cnm.estimates.A[:, cnm.estimates.idx_components[m2]].toarray()
A_TP = A_TP.reshape(dims + (-1,), order='F').transpose(2,0,1)
A_FN = estimates_gt.A[:, nm1].toarray()
A_FN = A_FN.reshape(dims + (-1,), order='F').transpose(2,0,1)
A_FP = A_onacid[:,cnm.estimates.idx_components[nm2]]
A_FP = A_FP.reshape(dims + (-1,), order='F').transpose(2,0,1)
plt.figure(figsize=(15, 12))
plt.imshow(cn_filter, vmin=lp, vmax=hp, cmap='viridis')
plt.colorbar();
for aa in A_TP:
plt.contour(aa, [0.05], colors='k');
for aa in A_FN:
plt.contour(aa, [0.05], colors='r');
for aa in A_FP:
plt.contour(aa, [0.25], colors='w');
cl = ['k', 'r', 'w']
lb = ['both', 'CNMF-E only', 'ring CNN only']
day = [mlines.Line2D([], [], color=cl[i], label=lb[i]) for i in range(3)]
plt.legend(handles=day, loc=3)
plt.axis('off');
plt.margins(0, 0);
plt.savefig(os.path.join(fld_name, 'ring_CNN_contours_gSig_3.pdf'), bbox_inches='tight', pad_inches=0)
In [ ]:
A_rej = cnm.estimates.A[:, cnm.estimates.idx_components_bad].toarray()
A_rej = A_rej.reshape(dims + (-1,), order='F').transpose(2,0,1)
plt.figure(figsize=(15, 15))
plt.imshow(cn_filter, vmin=lp, vmax=hp, cmap='viridis')
plt.title('Rejected Components')
for aa in A_rej:
plt.contour(aa, [0.05], colors='w');
In [ ]:
from caiman.utils.nn_models import create_LN_model
model_LN = create_LN_model(images, shape=opts.data['dims'] + (1,), n_channels=opts.ring_CNN['n_channels'],
width=opts.ring_CNN['width'], use_bias=opts.ring_CNN['use_bias'], gSig=gSig[0],
use_add=opts.ring_CNN['use_add'])
model_LN.load_weights(cnm.params.ring_CNN['path_to_model'])
In [ ]:
W = model_LN.get_weights()
In [ ]:
plt.figure(figsize=(10, 10))
plt.subplot(2,2,1); plt.imshow(np.squeeze(W[0][:,:,:,0])); plt.colorbar(); plt.title('Ring Kernel 1')
plt.subplot(2,2,2); plt.imshow(np.squeeze(W[0][:,:,:,1])); plt.colorbar(); plt.title('Ring Kernel 2')
plt.subplot(2,2,3); plt.imshow(np.squeeze(W[-1][:,:,0])); plt.colorbar(); plt.title('Multiplicative Layer 1')
plt.subplot(2,2,4); plt.imshow(np.squeeze(W[-1][:,:,1])); plt.colorbar(); plt.title('Multiplicative Layer 2');
In [ ]:
m1 = cm.load(fnames, var_name_hdf5='Y') # original data
m2 = cm.load(fname_mmap) # background subtracted data
m3 = m1 - m2 # estimated background
m4 = cm.movie(cnm.estimates.A[:,cnm.estimates.idx_components].dot(cnm.estimates.C[cnm.estimates.idx_components])).reshape(dims + (T,)).transpose(2,0,1)
# estimated components
nn = 0.01
mm = 1 - nn/4 # normalize movies by quantiles
m1 = (m1 - np.quantile(m1[:1000], nn))/(np.quantile(m1[:1000], mm) - np.quantile(m1[:1000], nn))
m2 = (m2 - np.quantile(m2[:1000], nn))/(np.quantile(m2[:1000], mm) - np.quantile(m2[:1000], nn))
m3 = (m3 - np.quantile(m3[:1000], nn))/(np.quantile(m3[:1000], mm) - np.quantile(m3[:1000], nn))
m4 = (m4 - np.quantile(m4[:1000], nn))/(np.quantile(m4[:1000], mm) - np.quantile(m4[:1000], nn))
m = cm.concatenate((cm.concatenate((m1.transpose(0,2,1), m3.transpose(0,2,1)), axis=2),
cm.concatenate((m2.transpose(0,2,1), m4), axis=2)), axis=1)
m[:3000].play(magnification=2, q_min=1, plot_text=True,
save_movie=True, movie_name=os.path.join(fld_name, 'movie.avi'))