get_ipython().magic('load_ext autoreload')
get_ipython().magic('autoreload 2')

import os
os.environ['MKL_NUM_THREADS'] = '1'
os.environ['OPENBLAS_NUM_THREADS'] = '1'

import bokeh.plotting as bpl
import cv2
import holoviews as hv
import logging
import matplotlib.pyplot as plt
from multiprocessing import Process, Queue
import numpy as np
from scipy.sparse import csr_matrix
from sys import platform
import tensorflow as tf
from time import time, sleep

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
from caiman.source_extraction.cnmf.online_cnmf import demix1p
from caiman.utils.nn_models import (fit_NL_model, create_LN_model, quantile_loss, rate_scheduler)
from caiman.utils.utils import download_demo

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

if tf.__version__[0] == '1':
    tf.enable_eager_execution()

"""
Explanation: Pipeline for real-time processing of microendoscopic data with CaImAn
This demo presents 3 approaches for processing microendoscopic data in real time using CaImAn. 
1. Sufficiently long initialization phase to identify all ROIs  followed by tracking
2. Short initalization phase followed by online processing using OnACID-E 
3. Short initalization phase followed by online processing using Ring-CNN+OnACID
All approached include:
- Motion Correction using the NoRMCorre algorithm
- Source extraction using a variant of the CNMF algorithm
- Deconvolution using the OASIS algorithm
OnACID-E and Ring-CNN further include
- Detection of new components
- Refinement of neural spatial footprints (and background in the case of OnACID-E)
End of explanation
"""


def get_iterator(device=0, fr=None):
    """
    device: device number (int) or filename (string) for reading from camera or file respectively
    fr: frame rate
    """
    if isinstance(device, int):  # capture from camera
        def capture_iter(device=device, fr=fr):
            cap = cv2.VideoCapture(device)
            if fr is not None:  # set frame rate
                cap.set(cv2.CAP_PROP_FPS, fr)
            while True:
                yield cv2.cvtColor(cap.read()[1], cv2.COLOR_BGR2GRAY)
        iterator = capture_iter(device, fr)
    else:  # read frame by frame from file
        iterator = cm.base.movies.load_iter(device, var_name_hdf5='Y')
    return iterator 

"""
Explanation: Below, we mention the code that can be used to read directly from the camera.<br/>
In this demo we simulate it by reading frame by frame from a file instead, so that everybody can execute the demo.
The download_demo function will download the file (if not already present) and store it inside your caiman_data/example_movies folder. We create an iterator over the file that returns the next imaged frame when calling next(iterator).<br/>
Here we define a function that creates and returns such an iterator. It works for both, a file or an imaging device.
End of explanation
"""


init_batch = 3000  # number of frames to use for initialization
T = 6000           # total number of frames 
fr = 10            # frame rate (Hz)

iterator = get_iterator(download_demo('blood_vessel_10Hz.mat'))

m = cm.movie(np.array([next(iterator) for t in range(init_batch)], dtype='float32'))

"""
Explanation: 1. Sufficiently long initialization phase to identify all ROIs  followed by tracking
Record for few minutes
The download_demo function will download the file (if not already present) and store it inside your caiman_data/example_movies folder.
End of explanation
"""


fname_init = m.save('init.mmap', order='C')

params_dict = {'fnames': fname_init,
               'fr': fr,
               'method_init': 'corr_pnr',
               'K': None,
               'gSig': (3, 3),
               'gSiz': (13, 13),
               'merge_thr': .65,
               'p': 1,
               'tsub': 1,
               'ssub': 1,
               'only_init': True,
               'nb': 0,
               'min_corr': .7,
               'min_pnr': 7,
               'normalize_init': False,
               'ring_size_factor': 1.4,
               'center_psf': True,
               'ssub_B': 2,
               'init_iter': 1,
               's_min': -10,
               'init_batch': init_batch,
               'init_method': 'cnmf',
               'use_dense': False,
               'motion_correct': True,     # flag for performing motion correction
               'gSig_filt': (3, 3)}        # size of high pass spatial filtering, used in 1p data             
opts = cnmf.params.CNMFParams(params_dict=params_dict)

cnm1 = cnmf.online_cnmf.OnACID(dview=None, params=opts)
cnm1.initialize_online(T=T)

"""
Explanation: Take a break from imaging to process recorded data
Taking a break to keep this demo simple. One could in parallel continue to save the otherwise "lost" frames to disk if one was not only intersted in the real-time experiment but post-analysis of the entire session
End of explanation
"""


# def init_from_other_pipeline(cnm, movie, A, C, T):
#     Yr = movie.reshape(len(m), -1).T
#     e = cnm.estimates
#     e.A = A
#     e.C = np.zeros_like(C)
#     e.bl = np.zeros(len(C))
#     e.c1 = np.zeros(len(C))
#     e.g = np.zeros((len(C), 1))
#     e.neurons_sn = np.zeros(len(C))
#     e.S = np.zeros_like(C)
#     e.lam = np.zeros(len(C))
#     for i,c in enumerate(C):
#         e.C[i], e.bl[i], e.c1[i], e.g[i], e.neurons_sn[i], e.S[i], e.lam[i] = (
#         cm.source_extraction.cnmf.deconvolution.constrained_foopsi(
#         c, p=1, bas_nonneg=False, noise_method='mean', fudge_factor=.97, optimize_g=5))
#     e.YrA = C-e.C
#     e.W, e.b = cm.source_extraction.cnmf.initialization.compute_W(
#         Yr, e.A, e.C, movie.shape[1:],
#         cnm.params.get('init', 'gSiz')[0] * cnm.params.get('init', 'ring_size_factor'),
#         ssub=cnm.params.get('init', 'ssub_B'))
#     cnm._prepare_object(Yr, T)
#     return cnm
#
# replace A,C with the results of your pipeline, and m with your initialization movie 
# A, C = cnm.estimates.Ab, cnm.estimates.noisyC[:cnm.N]
# cnm = cnmf.online_cnmf.OnACID(dview=None, params=opts)
# cnm = init_from_other_pipeline(cnm, m, A, C, T)

"""
Explanation: Alternatively we could initialize using any other pipeline that models the movie data as $Y=A C + B + N$ with spatial components $A$, temporal components $C$, background $B$, and noise $N$
End of explanation
"""


cnm1.t_read = []
cnm1.t_demix = []
cnm1.t_deconvolve = []
est = cnm1.estimates
ssub_B = cnm1.params.get('init', 'ssub_B') * cnm1.params.get('init', 'ssub')
for t in range(init_batch, T):
    # read frame
    t0 = time()
    frame = next(iterator)
    cnm1.t_read.append(time()-t0)
    # motion correct
    t0 = time()
    frame = cnm1.mc_next(t, frame.astype(np.float32))
    cnm1.t_motion.append(time()-t0)
    # get noisy fluorescence value via NNLS (project data on shapes & demix)
    t0 = time()
    C_in = est.noisyC[:cnm1.M, t - 1].copy()
    est.C_on[:cnm1.M, t], est.noisyC[:cnm1.M, t] = demix1p(
        frame.ravel(order='F'), est.Ab, C_in, est.AtA,
        Atb=est.Atb, AtW=est.AtW, AtWA=est.AtWA,
        iters=3, groups=est.groups, ssub_B=ssub_B, 
        downscale_matrix=est.downscale_matrix if ssub_B > 1 else None)
    cnm1.t_demix.append(time()-t0)
    # denoise & deconvolve
    t0 = time()
    if cnm1.params.get('preprocess', 'p'):
        for i, o in enumerate(est.OASISinstances):
            o.fit_next(est.noisyC[i, t])
            est.C_on[i, t - o.get_l_of_last_pool() +
                      1: t + 1] = o.get_c_of_last_pool()
    cnm1.t_deconvolve.append(time()-t0)
    # add code to display whatever you want in order to guide the closed-loop experiment below
    # e.g. print the indices of neurons that just spiked
#     print('\r', ' '*200, end="\r")
#     for i, o in enumerate(est.OASISinstances):
#         if o.get_l_of_last_pool() == 1:
#             print(i, end=' ')

del iterator

"""
Explanation: Start real-time processing
End of explanation
"""


# calculate time one would have to wait for next frame to arrive if reading directly from camera
t_all = np.cumsum(cnm1.t_read) + np.cumsum(cnm1.t_motion) + np.cumsum(cnm1.t_demix) + np.cumsum(cnm1.t_deconvolve)

t_wait=[1]
t_wait_total=0
for i in range(1, T-init_batch):
    t_wait.append(max(i/fr - t_all[i-1]-t_wait_total, 0))
    t_wait_total += t_wait[-1]
realtime = np.array(t_wait)>0
print('%g%s processed in real time. %g/%g frames' % 
      (100 * realtime.sum() / (T-init_batch), '%', realtime.sum(), T-init_batch))

plt.figure(figsize=(12,4))
for i, f in enumerate((lambda a: 1000*np.array(a), np.cumsum)):
    plt.subplot(1,2,1+i)
    plt.stackplot(np.arange(T-init_batch), f(cnm1.t_read), f(cnm1.t_motion),
                  f(cnm1.t_demix), f(cnm1.t_deconvolve))
    plt.gca().add_artist(plt.legend(labels=['read', 'motion', 'demix', 'deconvolve'], loc=2))
    plt.title('Processing time allocation')
    plt.xlabel('Frame #')
    plt.ylabel(('Processing time per frame [ms]', 'Cumulative processing time [s]')[i])
    if i==0:
        plt.fill_between(range(T-init_batch),[0]*(T-init_batch),
                 [plt.ylim()[1]]*(T-init_batch), where=realtime, 
                 color='y', alpha=.1, edgecolor='y', zorder=-11, label='real time')
        plt.gca().add_artist(plt.legend())

cnm1.estimates.A = cnm1.estimates.Ab
cnm1.estimates.C = cnm1.estimates.C_on
cnm1.estimates.YrA = cnm1.estimates.noisyC-cnm1.estimates.C

cn, pnr = cm.summary_images.correlation_pnr(cm.load(download_demo('blood_vessel_10Hz.mat'), var_name_hdf5='Y'), gSig=3, swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile
cnm1.estimates.coordinates = None
cnm1.estimates.plot_contours_nb(img=cn, thr=.6);

cnm1.estimates.nb_view_components(img=cn, denoised_color='red');

"""
Explanation: Plot results
End of explanation
"""


init_batch = 500  # number of frames to use for initialization
T = 6000          # total number of frames 
fr = 10           # frame rate (Hz)

iterator = get_iterator(download_demo('blood_vessel_10Hz.mat'))

m = cm.movie(np.array([next(iterator) for t in range(init_batch)], dtype='float32'))

"""
Explanation: 2. Short initalization phase followed by online processing using OnACID-E
Record for some seconds
The download_demo function will download the file (if not already present) and store it inside your caiman_data/example_movies folder.
End of explanation
"""


fname_init = m.save('init.mmap', order='C')

params_dict = {'fnames': fname_init,
               'fr': fr,
               'method_init': 'corr_pnr',
               'K': None,
               'gSig': (3, 3),
               'gSiz': (13, 13),
               'merge_thr': .65,
               'p': 1,
               'tsub': 1,
               'ssub': 1,
               'only_init': True,
               'nb': 0,
               'min_corr': .65,
               'min_pnr': 6,
               'normalize_init': False,
               'ring_size_factor': 1.4,
               'center_psf': True,
               'ssub_B': 2,
               'init_iter': 1,
               's_min': -10,
               'init_batch': init_batch,
               'init_method': 'cnmf',
               'batch_update_suff_stat': True,
               'update_freq': 200,
               'expected_comps': 600,
               'update_num_comps': True,
               'rval_thr': .55,
               'thresh_fitness_raw': -130,
               'thresh_fitness_delta': -20,
               'min_SNR': 2.5,
               'min_num_trial': 5,
               'max_num_added': 5,
               'use_corr_img': False,
               'use_dense': False,
               'motion_correct': True,     # flag for performing motion correction
               'gSig_filt': (3, 3),        # size of high pass spatial filtering, used in 1p data
               'use_cnn': False}
opts = cnmf.params.CNMFParams(params_dict=params_dict)

cnm2 = cnmf.online_cnmf.OnACID(dview=None, params=opts)
cnm2.initialize_online(T=T)

"""
Explanation: Take a break from imaging to process recorded data
Taking a break to keep this demo simple. One could in parallel continue to save the otherwise "lost" frames to disk if one was not only intersted in the real-time experiment but post-analysis of the entire session
End of explanation
"""


if platform in ('linux', 'darwin'):
    """ This code assumes you are running on a Unix-based computing system (like Linux or macOS)
        that allows to start a Process uing a fork() call, NOT Windows. For the latter see the else call below
        https://docs.python.org/3.7/library/multiprocessing.html#the-spawn-and-forkserver-start-methods"""
    q = Queue()
    cnm2.t_wait = []
    cnm2.t_motion = []
    cnm2.t_fit = []
    realtime = []

    def append_to_queue(q, init_batch, T, fr):
        t_start = time()
        for t in range(init_batch, T):
            # read frame and append to queue
            frame = next(iterator)
            q.put(frame)
            sleep(max(0, (t+1-init_batch)/fr - time()+t_start))

    producer = Process(target=append_to_queue, args=(q,init_batch,T,fr))
    producer.start()

    # process first frame
    t = init_batch
    # read form queue (wait for next frame if empty)
    t0=time()
    frame = q.get()
    cnm2.t_wait.append(time()-t0)
    t_start = time() + (1-init_batch)/fr
    # motion correct
    t0 = time()
    frame = cnm2.mc_next(t, frame)
    cnm2.t_motion.append(time()-t0)
    # fit
    t0 = time()
    cnm2.fit_next(t, frame.ravel(order='F'))
    cnm2.t_fit.append(time()-t0)
    rt = time() <= t_start + t/fr
    realtime.append(rt)

    # process remaining frames
    for t in range(init_batch+1, T):
        # read form queue (wait for next frame if empty)
        t0=time()
        frame = q.get()
        cnm2.t_wait.append(time()-t0)
        # motion correct
        t0 = time()
        frame = cnm2.mc_next(t, frame)
        cnm2.t_motion.append(time()-t0)
        # fit
        t0 = time()
        cnm2.fit_next(t, frame.ravel(order='F'))
        cnm2.t_fit.append(time()-t0)
        rt = time() <= t_start + t/fr
        realtime.append(rt)
        # add code to display whatever you want in order to guide the closed-loop experiment below
        print('Realtime: ' + ("\x1b[32mTrue\x1b[0m" if rt else "\x1b[31mFalse\x1b[0m"), end="  \r", flush=True)
else:

    """Windows requires as workaround to move some code to the external file 'win.py'.
       The new Process is being spawned not forked. Expect performance to be somewhat impeded."""
    print('Windows lacks os.fork(). Expect performance to be somewhat impeded.')
    import win

    if __name__ == '__main__':
        q = Queue()
        cnm2.t_wait = []
        cnm2.t_motion = []
        cnm2.t_fit = []
        realtime = []
        producer = Process(target=win.append_to_queue, args=(q,init_batch,T,fr))
        producer.start()

        # process first frame
        t = init_batch
        # read form queue (wait for next frame if empty)
        t0=time()
        frame = q.get()
        cnm2.t_wait.append(time()-t0)
        t_start = time() + (1-init_batch)/fr
        # motion correct
        t0 = time()
        frame = cnm2.mc_next(t, frame)
        cnm2.t_motion.append(time()-t0)
        # fit
        t0 = time()
        cnm2.fit_next(t, frame.ravel(order='F'))
        cnm2.t_fit.append(time()-t0)
        rt = time() <= t_start + t/fr
        realtime.append(rt)

        #process remaining frames
        for t in range(init_batch+1, T):
            # read form queue (wait for next frame if empty)
            t0=time()
            frame = q.get()
            cnm2.t_wait.append(time()-t0)
            # motion correct
            t0 = time()
            frame = cnm2.mc_next(t, frame)
            cnm2.t_motion.append(time()-t0)
            # fit
            t0 = time()
            cnm2.fit_next(t, frame.ravel(order='F'))
            cnm2.t_fit.append(time()-t0)
            rt = time() <= t_start + t/fr
            realtime.append(rt)
            # add code to display whatever you want in order to guide the closed-loop experiment below
            print('Realtime: ' + ("\x1b[32mTrue\x1b[0m" if rt else "\x1b[31mFalse\x1b[0m"), end="  \r", flush=True)
        producer.join()

del iterator

"""
Explanation: Start real-time processing
The spatial footprints and background parameters are updated every 200 frames. This results in a long processing time for the individual frames for which these updates occur. Thus instead of the usual waiting for the camera to provide the next image, there are few frames for which the images are acquired faster than processed. We use a separate Process to acquire and add the images to a FIFO Queue at regular time intervals. The main process reads the next image from this Queue and waits for it if the Queue is empty.
End of explanation
"""


print('%g%s processed in real time. %g/%g frames' % 
      (100 * np.sum(realtime) / (T-init_batch), '%', np.sum(realtime), T-init_batch))

plt.figure(figsize=(12,4))
cnm2.t_read = np.zeros(T-init_batch)
for i, f in enumerate((lambda a: 1000*np.array(a), np.cumsum)):
    plt.subplot(1,2,1+i)
    plt.stackplot(np.arange(len(cnm2.t_fit)), f(cnm2.t_read), f(cnm2.t_motion),
                  f(np.array(cnm2.t_fit) - np.array([cnm2.t_detect, cnm2.t_shapes, cnm2.t_stat]).sum(0)),
                  f(cnm2.t_detect), f(cnm2.t_shapes)+f(cnm2.t_stat))
    plt.gca().add_artist(plt.legend(labels=['read', 'motion', 'process', 'detect', 'shapes'], loc=2))
    plt.title('Processing time allocation')
    plt.xlabel('Frame #')
    plt.ylabel(('Processing time per frame [ms]', 'Cumulative processing time [s]')[i])
    if i==0:
        plt.ylim(0, 100)
        plt.fill_between(range(T-init_batch),[0]*(T-init_batch),
                 [100]*(T-init_batch), where=realtime, 
                 color='y', alpha=.1, edgecolor='y', zorder=-11, label='real time')
        plt.gca().add_artist(plt.legend())

cnm2.estimates.A = cnm2.estimates.Ab
cnm2.estimates.C = cnm2.estimates.C_on[:cnm2.N]
cnm2.estimates.YrA = cnm2.estimates.noisyC[:cnm2.N]-cnm2.estimates.C

cn, pnr = cm.summary_images.correlation_pnr(cm.load(download_demo('blood_vessel_10Hz.mat'), var_name_hdf5='Y'), gSig=3, swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile
cnm2.estimates.coordinates = None
cnm2.estimates.plot_contours_nb(img=cn, thr=.6);

cnm2.estimates.nb_view_components(img=cn, denoised_color='red');

"""
Explanation: Plot results
End of explanation
"""


init_batch = 500  # number of frames to use for initialization
T = 6000          # total number of frames 
fr = 10           # frame rate (Hz)

iterator = get_iterator(download_demo('blood_vessel_10Hz.mat'))

m = cm.movie(np.array([next(iterator) for t in range(init_batch)], dtype='float32'))

"""
Explanation: 3. Short initalization phase followed by online processing using Ring-CNN+OnACID
Record for some seconds
The download_demo function will download the file (if not already present) and store it inside your caiman_data/example_movies folder.
End of explanation
"""


fname_init = m.save('init.mmap', order='C')

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

params_dict = {'fnames': fname_init,
               'var_name_hdf5': 'Y',                                # name of variable inside mat file where the data is stored
               'fr': fr,                                            # frame rate (Hz)
               'decay_time': 0.5,                                   # approximate length of transient event in seconds
               'gSig': gSig,
               'p': 1,                                              # 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': True,                              # 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': 700,                                # 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': 1,                                          # number of total passes over the data
               'stop_detection': False,                              # 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)

"""
Explanation: Take a break from imaging to process recorded data
Taking a break to keep this demo simple. One could in parallel continue to save the otherwise "lost" frames to disk if one was not only intersted in the real-time experiment but post-analysis of the entire session
End of explanation
"""


cnm3 = cnmf.online_cnmf.OnACID(params=opts)

if cnm3.params.get('ring_CNN', 'loss_fn') == 'pct':
    loss_fn = quantile_loss(cnm3.params.get('ring_CNN', 'pct'))
else:
    loss_fn = cnm3.params.get('ring_CNN', 'loss_fn')
if cnm3.params.get('ring_CNN', 'lr_scheduler') is None:
    sch = None
else:
    sch = rate_scheduler(*cnm3.params.get('ring_CNN', 'lr_scheduler'))
model_LN = create_LN_model(m, shape=m.shape[1:] + (1,),
                           n_channels=cnm3.params.get('ring_CNN', 'n_channels'),
                           lr=cnm3.params.get('ring_CNN', 'lr'),
                           gSig=cnm3.params.get('init', 'gSig')[0],
                           loss=loss_fn, width=cnm3.params.get('ring_CNN', 'width'),
                           use_add=cnm3.params.get('ring_CNN', 'use_add'),
                           use_bias=cnm3.params.get('ring_CNN', 'use_bias'))
if cnm3.params.get('ring_CNN', 'reuse_model'):
    model_LN.load_weights(cnm3.params.get('ring_CNN', 'path_to_model'))
else:
    model_LN, history, path_to_model = fit_NL_model(
        model_LN, m, epochs=cnm3.params.get('ring_CNN', 'max_epochs'),
        patience=cnm3.params.get('ring_CNN', 'patience'), schedule=sch)
    cnm3.params.set('ring_CNN', {'path_to_model': path_to_model})

"""
Explanation: Train or load Ring CNN
End of explanation
"""


cnm3.initialize_online(T=T, model_LN=model_LN)

# if no GPU is available prediction is faster when the CNN is converted into a sparse matrix 
if tf.test.is_gpu_available():
    predict = model_LN.predict
    # dummy prediction to initialize CNN-ring model model_LN
    predict(m[:1, ..., None])
else:
    dims = m[0].shape
    kernel_size = model_LN.layers[1].weights[0].shape[0]
    kernel_half = (kernel_size-1)//2
    w_1 = model_LN.layers[1].weights[0].numpy().squeeze().transpose(2,0,1)
    w_1 = w_1.reshape(len(w_1), -1)
    w_2 = model_LN.layers[2].weights[0].numpy()
    data = []
    indices = []
    indptr = [0]
    for i in range(dims[0]):
        for j in range(dims[1]):
            tmp = np.zeros((dims[0]+kernel_size, dims[1]+kernel_size), dtype=np.float32)
            tmp[i:i+kernel_size, j:j+kernel_size] = w_2[i, j].dot(w_1).reshape(kernel_size, kernel_size)
            tmp = tmp[kernel_half:kernel_half+dims[0], kernel_half:kernel_half+dims[1]].ravel()
            ind = np.where(tmp!=0)[0]
            newdata = list(tmp[ind])
            data += newdata
            indices += list(ind)
            indptr += [len(newdata) + indptr[-1]]
    w = csr_matrix((data, indices, indptr)).todia()

    def predict(frame):
        return w.dot(frame.ravel()).reshape(dims)

"""
Explanation: Initialize OnACID
End of explanation
"""


cnm3.t_read = []
cnm3.t_bkgrd = []
cnm3.t_motion = []
cnm3.t_fit = []
for t in range(init_batch, T):
    # read frame
    t0 = time()
    frame = next(iterator)
    cnm3.t_read.append(time()-t0)
    # remove background
    t0 = time()
    frame = np.maximum(frame - np.squeeze(predict(
        frame.astype(np.float32)[None,...,None])), 0)
    cnm3.t_bkgrd.append(time()-t0)
    # motion correct
    t0 = time()
    frame = cnm3.mc_next(t, frame.astype(np.float32))
    cnm3.t_motion.append(time()-t0)
    # fit
    t0 = time()
    cnm3.fit_next(t, frame.ravel(order='F'))
    cnm3.t_fit.append(time()-t0)
    # add code to display whatever you want in order to guide the closed-loop experiment below

del iterator

"""
Explanation: Start real-time processing
End of explanation
"""


# calculate time one would have to wait for next frame to arrive if reading directly from camera
t_all = np.cumsum(cnm3.t_read) + np.cumsum(cnm3.t_bkgrd) + np.cumsum(cnm3.t_motion) + np.cumsum(cnm3.t_fit)

t_wait=[1]
t_wait_total=0
for i in range(1, T-init_batch):
    t_wait.append(max(i/fr - t_all[i-1]-t_wait_total, 0))
    t_wait_total += t_wait[-1]
realtime = np.array(t_wait)>0
print('%g%s processed in real time. %g/%g frames' % 
      (100 * realtime.sum() / (T-init_batch), '%', realtime.sum(), T-init_batch))

plt.figure(figsize=(12,4))
for i, f in enumerate((lambda a: 1000*np.array(a), np.cumsum)):
    plt.subplot(1,2,1+i)
    plt.stackplot(np.arange(len(cnm3.t_fit)), f(cnm3.t_read), f(cnm3.t_bkgrd), f(cnm3.t_motion),
                  f(np.array(cnm3.t_fit) - np.array([cnm3.t_detect, cnm3.t_shapes, cnm3.t_stat]).sum(0)),
                  f(cnm3.t_detect), f(cnm3.t_shapes)+f(cnm3.t_stat))
    plt.gca().add_artist(plt.legend(labels=['read', 'background', 'motion', 'process', 'detect', 'shapes'], loc=2))
    plt.title('Processing time allocation')
    plt.xlabel('Frame #')
    plt.ylabel(('Processing time per frame [ms]', 'Cumulative processing time [s]')[i])
    if i==0:
        plt.fill_between(range(T-init_batch),[0]*(T-init_batch),
                 [plt.ylim()[1]]*(T-init_batch), where=realtime, 
                 color='y', alpha=.1, edgecolor='y', zorder=-11, label='real time')
        plt.gca().add_artist(plt.legend())

cnm3.estimates.A = cnm3.estimates.Ab[:,gnb:]
cnm3.estimates.C = cnm3.estimates.C_on[gnb:cnm3.M]
cnm3.estimates.YrA = cnm3.estimates.noisyC[gnb:cnm3.M]-cnm3.estimates.C

cn, pnr = cm.summary_images.correlation_pnr(cm.load(download_demo('blood_vessel_10Hz.mat'), var_name_hdf5='Y'), gSig=3, swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile
cnm3.estimates.coordinates = None
cnm3.estimates.plot_contours_nb(img=cn, thr=.6);

cnm3.estimates.nb_view_components(img=cn,denoised_color='red');

"""
Explanation: Plot results
End of explanation
"""


