Online CNMF-E

This demo shows an example of doing online analysis on one-photon data. We compare offline and online approaches. The dataset used is courtesy of the Miniscope project.



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 scipy

logging.basicConfig(format=
                          "%(relativeCreated)12d [%(filename)s:%(funcName)10s():%(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.motion_correction import MotionCorrect
from caiman.utils.utils import download_demo
import matplotlib.pyplot as plt
from caiman.utils.visualization import nb_inspect_correlation_pnr
import holoviews as hv
import bokeh.plotting as bpl

bpl.output_notebook()
hv.notebook_extension('bokeh')

Select file(s) to be processed

The download_demo function will download the specific file for you and return the complete path to the file which will be stored in your caiman_data directory. If you adapt this demo for your data make sure to pass the complete path to your file(s). Remember to pass the fnames variable as a list. Note that the memory requirement of the offline CNMF-E algorithm are much higher compared to the standard CNMF algorithm. One of the benefits of the online approach is the reduced memory requirements.



In [ ]:

    
fnames = [download_demo('msCam13.avi')]

Batch (offline) approach

We start with motion correction and then proceed with the source extraction using the CNMF-E algorithm. For a detailed 1p demo check demo_pipeline_cnmfE.ipynb.



In [ ]:

    
# motion correction parameters
motion_correct = True            # flag for performing motion correction
pw_rigid = False                 # flag for performing piecewise-rigid motion correction (otherwise just rigid)
gSig_filt = (7, 7)               # size of high pass spatial filtering, used in 1p data
max_shifts = (20, 20)            # maximum allowed rigid shift
border_nan = 'copy'              # replicate values along the boundaries

mc_dict = {
    'pw_rigid': pw_rigid,
    'max_shifts': max_shifts,
    'gSig_filt': gSig_filt,
    'border_nan': border_nan
}

opts = cnmf.params.CNMFParams(params_dict=mc_dict)



In [ ]:

    
#%% start a cluster for parallel processing (if a cluster already exists it will be closed and a new session will be opened)
if 'dview' in locals():
    cm.stop_server(dview=dview)
c, dview, n_processes = cm.cluster.setup_cluster(
    backend='local', n_processes=None, single_thread=False)



In [ ]:

    
mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion'))
mc.motion_correct(save_movie=True)

inspect motion correction results



In [ ]:

    
inspect_results = False
if inspect_results:
    cm.concatenate((cm.load(fnames), cm.load(mc.mmap_file)), axis=1).play()
#plt.figure(); plt.plot(mc.shifts_rig); plt.legend(['x-shifts', 'y-shifts'])

The motion correction results look good. We then proceed with memory mapping and checking the correlation/pnr images.



In [ ]:

    
from time import time
fname_new = cm.save_memmap(mc.mmap_file, base_name='memmap_', order='C',
                           border_to_0=0, dview=dview)
Yr, dims, T = cm.load_memmap(fname_new)
images = Yr.T.reshape((T,) + dims, order='F')

Inspect correlation and PNR images to set relevant thresholds



In [ ]:

    
gSig = (6, 6)

cn_filter, pnr = cm.summary_images.correlation_pnr(images[::max(T//1000, 1)], gSig=gSig[0], swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile
# inspect the summary images and set the parameters

nb_inspect_correlation_pnr(cn_filter, pnr)

Set parameters for source extraction

From the images above we select min_pnr = 10 and min_corr = 0.8. We pass these alongside the other parameters needed for offline 1p processing.



In [ ]:

    
min_pnr = 10
min_corr = 0.8
rf = 48                                        # half size of each patch
stride = 8                                     # amount of overlap between patches     
ssub = 1                                       # spatial downsampling factor   
decay_time = 0.4                               # length of typical transient (in seconds)
fr = 10                                        # imaging rate (Hz) 
gSig = (6, 6)                                  # expected half size of neurons
gSiz = (15, 15)                                # half size for neuron bounding box   
p = 0                                          # order of AR indicator dynamics
min_SNR = 1.5                                  # minimum SNR for accepting new components
rval_thr = 0.85                                # correlation threshold for new component inclusion
merge_thr = 0.65                               # merging threshold
K = None                                       # initial number of components

cnmfe_dict = {'fnames': fnames,
              'fr': fr,
              'decay_time': decay_time,
              'method_init': 'corr_pnr',
              'gSig': gSig,
              'gSiz': gSiz,
              'rf': rf,
              'stride': stride,
              'p': p,
              'nb': 0,
              'ssub': ssub,
              'min_SNR': min_SNR,
              'min_pnr': min_pnr,
              'min_corr': min_corr,
              'bas_nonneg': False,
              'center_psf': True,
              'rval_thr': rval_thr,
              'only_init': True,
              'merge_thr': merge_thr,
              'K': K}
opts.change_params(cnmfe_dict);



In [ ]:

    
from time import time
t1 = -time()
cnm = cnmf.CNMF(n_processes=n_processes, dview=dview, params=opts)
cnm.fit(images)
t1 += time()

View the results



In [ ]:

    
cnm.estimates.plot_contours_nb(img=pnr)



In [ ]:

    
cnm.estimates.hv_view_components(img=cn_filter)

Show a movie with the results



In [ ]:

    
cnm.estimates.play_movie(images, magnification=0.75, include_bck=False)

Online Processing

Now try the online approach. The idea behind the online algorithm is simple:

First initialize the estimates by running the batch (offline) algorithm in small subset.
Then process each frame as it arrives. The processing consists of:
- Motion correct the new frame
- Extract the activity of existing neurons at this frame, and neuropil
- Search for new neurons that appear in this frame and have not been detected earlier.
Periodically update shapes of existing neurons and background model.

Setup additional parameters for online processing



In [ ]:

    
from copy import deepcopy
online_opts = deepcopy(cnm.params)



In [ ]:

    
rf = 48                                        # half size of patch (used only during initialization)
stride = 8                                     # overlap between patches (used only during initialization) 
ssub = 1                                       # spatial downsampling factor (during initialization)
ds_factor = 2*ssub                             # spatial downsampling factor (during online processing)         
ssub_B = 4                                     # background downsampling factor (use that for faster processing)
gSig = (10//ds_factor, 10//ds_factor)          # expected half size of neurons
gSiz = (22//ds_factor, 22//ds_factor)
sniper_mode = False                            # flag using a CNN to detect new neurons (o/w space correlation is used)
init_batch = 300                               # number of frames for initialization (presumably from the first file)
expected_comps = 500                           # maximum number of expected components used for memory pre-allocation (exaggerate here)
dist_shape_update = False                      # flag for updating shapes in a distributed way
min_num_trial = 5                              # number of candidate components per frame     
K = None                                       # initial number of components
epochs = 2                                     # number of passes over the data
show_movie = False                             # show the movie with the results as the data gets processed
use_corr_img = True                            # flag for using the corr*pnr image when searching for new neurons (otherwise residual)

online_dict = {'epochs': epochs,
               'nb': 0,
               'ssub': ssub,
               'ssub_B': ssub_B,
               'ds_factor': ds_factor,                                   # ds_factor >= ssub should hold
               'gSig': gSig,
               'gSiz': gSiz,
               'gSig_filt': (3, 3),
               'min_corr': min_corr,
               'bas_nonneg': False,
               'center_psf': True,
               'max_shifts_online': 20,
               'rval_thr': rval_thr,
               'motion_correct': True,
               'init_batch': init_batch,
               'only_init': True,
               'init_method': 'cnmf',
               'normalize_init': False,
               'update_freq': 200,
               'expected_comps': expected_comps,
               'sniper_mode': sniper_mode,                               # set to False for 1p data       
               'dist_shape_update' : dist_shape_update,
               'min_num_trial': min_num_trial,
               'epochs': epochs,
               'use_corr_img': use_corr_img,
               'show_movie': show_movie}
online_opts.change_params(online_dict);



In [ ]:

    
cnm_online = cnmf.online_cnmf.OnACID(params=online_opts, dview=dview)
cnm_online.fit_online()



In [ ]:

    
#images = cm.load(fnames[0], subindices=slice(0,1000))
#Cn, pnr = cm.summary_images.correlation_pnr(images[::1], gSig=gSig[0], swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile
cnm_online.estimates.nb_view_components(img=pnr, denoised_color='red');



In [ ]:

    
cnm_online.estimates.plot_contours_nb(img=pnr)

Plot timing

The plot below shows the time spent on each part of the algorithm (motion correction, tracking of current components, detect new components, update shapes) for each frame. Note that if you displayed a movie while processing the data (show_movie=True) the time required to generate this movie will be included here.



In [ ]:

    
show_cumulative = True
#if show_cumulative:
T_init = np.array([cnm_online.t_init] + [0]*(epochs*T-1))
T_motion = 1e3*np.array([0]*init_batch + cnm_online.t_motion)/1e3
T_detect = 1e3*np.array([0]*init_batch + cnm_online.t_detect)/1e3
T_shapes = 1e3*np.array([0]*init_batch + cnm_online.t_shapes)/1e3
T_online = 1e3*np.array([0]*init_batch + cnm_online.t_online)/1e3 - T_motion - T_detect - T_shapes
plt.figure()
plt.stackplot(np.arange(len(T_motion)), np.cumsum(T_init), np.cumsum(T_motion), np.cumsum(T_online), np.cumsum(T_detect), np.cumsum(T_shapes))
plt.legend(labels=['init', 'motion', 'process', 'detect', 'shapes'], loc=2)
for i in range(epochs - 1):
    plt.plot([(i+1)*T, (i+1)*T], [0, np.array(cnm_online.t_online).sum()+cnm_online.t_init], '--k')
plt.title('Processing time allocation')
plt.xlabel('Frame #')
plt.ylabel('Processing time [ms]')
#plt.ylim([0, 1.2e3*np.percentile(np.array(cnm_online.t_online), 90)]);



In [ ]:

    
cnm_online.estimates.play_movie(imgs=images, magnification=0.75, include_bck=False)

Clean up and compare two approaches

Even though the online algorithm screens any new components, we can still perform the quality tests to filter out any false positive components. To do that, we first need to apply the inferred shifts to the original data in order to have the whole registered dataset in memory mapped form.



In [ ]:

    
if online_opts.online['motion_correct']:
    shifts = cnm_online.estimates.shifts[-cnm_online.estimates.C.shape[-1]:]
    if not opts.motion['pw_rigid']:
        memmap_file = cm.motion_correction.apply_shift_online(images, shifts,
                                                              save_base_name='MC')
    else:
        mc = MotionCorrect(fnames, dview=dview, **online_opts.get_group('motion'))

        mc.y_shifts_els = [[sx[0] for sx in sh] for sh in shifts]
        mc.x_shifts_els = [[sx[1] for sx in sh] for sh in shifts]
        memmap_file = mc.apply_shifts_movie(fnames, rigid_shifts=False,
                                            save_memmap=True,
                                            save_base_name='MC')
else:  # To do: apply non-rigid shifts on the fly
    memmap_file = images.save(fnames[0][:-4] + 'mmap')
cnm_online.mmap_file = memmap_file
Yr_online, dims, T = cm.load_memmap(memmap_file)

#cnm_online.estimates.dview=dview
#cnm_online.estimates.compute_residuals(Yr=Yr_online)
images_online = np.reshape(Yr_online.T, [T] + list(dims), order='F')
min_SNR = 2  # peak SNR for accepted components (if above this, acept)
rval_thr = 0.85  # space correlation threshold (if above this, accept)
use_cnn = False  # use the CNN classifier
cnm_online.params.change_params({'min_SNR': min_SNR,
                                'rval_thr': rval_thr,
                                'use_cnn': use_cnn})

cnm_online.estimates.evaluate_components(images_online, cnm_online.params, dview=dview)
cnm_online.estimates.Cn = pnr



In [ ]:

    
cnm_online.estimates.plot_contours_nb(img=pnr, idx=cnm_online.estimates.idx_components)



In [ ]:

    
cnm_online.estimates.hv_view_components(img=pnr, idx=cnm_online.estimates.idx_components,
                                        denoised_color='red')



In [ ]:

    
cnm_online.estimates.hv_view_components(img=pnr, idx=cnm_online.estimates.idx_components_bad,
                                        denoised_color='red')

Difference in inferred shifts

Accurate motion correction is important for the online algorithm. Below we plot the difference in the estimated shifts between the two approaches. Note that the online shifts have been rescaled by a factor of ds_factor.



In [ ]:

    
plt.plot(np.array(mc.shifts_rig) - ds_factor*np.array(cnm_online.estimates.shifts[:1000]));
plt.legend(['x-shifts', 'y-shifts']);
plt.title('Difference between offline and online shifts')
plt.xlabel('Frame #')
plt.ylabel('pixels')

Constant shifts in the FOV will not significantly affect the results. What is most important is deviatons.



In [ ]:

    
np.std(np.array(mc.shifts_rig) - ds_factor*np.array(cnm_online.estimates.shifts[:1000]), axis=0)

The standard deviation is at a subpixel level (although it can still be significant). The high degree of similarity can also be seen from the correlation between the shifts of the two approaches.



In [ ]:

    
np.corrcoef(np.array(mc.shifts_rig).T, np.array(cnm_online.estimates.shifts[:1000]).T)