<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Analyzing dendritic data with CaImAn
This notebook shows an example on how to analyze two-photon dendritic data with CaImAn. It follows closely the other notebooks.
Step1:  Selecting the data
We provide an example file from mouse barrel cortex, courtesy of Clay Lacefield and Randy Bruno, Columbia University. The download_demo command will automatically download the file and store it in your caiman_data folder the first time you run it. To use the demo in your own dataset you can set
Step2:  Viewing the data
Prior to analysis you can view the data by setting display_movie = False. The example file provided here is small, with FOV size 128 x 128. To view it larger we set magnification = 4 in the play command. This will result to a movie of resolution 512 x 512. For larger files, you will want to reduce this magnification factor (for memory and display purposes).
Step3:  Parameter setting
First we'll perform motion correction followed by CNMF.
One of the main differences when analyzing dendritic data, is that dendritic segments are not localized and can traverse significant parts of the FOV. They can also be spatially non-contiguous. To capture this property, CaImAn has two methods for initializing CNMF
Step4:  Motion correction
First we create a motion correction object and then perform the registration. The provided dataset has already been registered with a different method so we do not expect to see a lot of motion.
Step5:  Display registered movie next to the original
Step6:  Memory mapping
Step7:  Now run CNMF with the chosen initialization method
Note that in the parameter setting we have set only_init = True so only the initialization will be run.
Step8:  Display components
And compute the correlation image as a reference image
Step9:  Plot a montage array of all the components
We can also plot a montage of all the identified spatial footprints.
Step10:  Run CNMF to sparsify the components
The components found by the initialization have captured a lot of structure but they are not sparse. Using the refit command will sparsify them.
Step11:  Now plot the components again
Step12:  The resulting components are significantly more sparse and capture a lot of the data structure.
Play a movie with the results
The movie will show three panels with the original, denoised and residual movies respectively. 
The flag include_bck = False will remove the estimated background from the original and denoised panels.
To include it set it to True.
Step13:  Save the object
You can save the results of the analysis for future use and/or to load them on the GUI.
Step14:  Now try CNMF with patches
As mentioned above, the file shown here is small enough to be processed at once. For larger files patches and/or downsampling can be used. We show how to use patches below. Note that unlike in somatic imaging, for dendritic imaging we want the patches to be large and have significant overlap so enough structure of the dendritic segments is captured in the overlapping region. This will help with merging the components. However, a too large patch size can result to memory issues due to large amount of data loaded concurrently onto memory.
Step15:  Change some parameters to enable patches
To enable patches the rf parameter should be changed from None to some integer value. Moreover the amount of overlap controlled by stride_cnmf should also change. Finally, since the parameter K specifies components per patch it should be reduced compared to its value when operating on the whole FOV at once. It is important to experiment with different values for these parameters since results can sensitive, due to the non stereotypical shapes and activity patterns of dendrites.
Step16:  Refit to sparsify
As before we can pass the results through the CNMF refit function that can yield more sparse components.
Step17:  Compare the two approaches
To identify how similar components and the two approaches yield, we can compute the cross-correlation matrix between the spatial footprints derived from the two approaches.
Step18:  Stop the cluster when you're done
<ASSISTANT_TASK:>
Python Code:
import cv2
import glob
import logging
import matplotlib.pyplot as plt
import numpy as np
import os
try:
    cv2.setNumThreads(0)
except():
    pass
try:
    if __IPYTHON__:
        get_ipython().magic('load_ext autoreload')
        get_ipython().magic('autoreload 2')
except NameError:
    pass
import caiman as cm
from caiman.motion_correction import MotionCorrect
from caiman.source_extraction.cnmf import cnmf as cnmf
from caiman.source_extraction.cnmf import params as params
from caiman.utils.utils import download_demo
from caiman.utils.visualization import nb_view_patches
from skimage.util import montage
import bokeh.plotting as bpl
import holoviews as hv
bpl.output_notebook()
hv.notebook_extension('bokeh')
Explanation: Analyzing dendritic data with CaImAn
This notebook shows an example on how to analyze two-photon dendritic data with CaImAn. It follows closely the other notebooks.
End of explanation
fnames = [download_demo('data_dendritic.tif')]
Explanation: Selecting the data
We provide an example file from mouse barrel cortex, courtesy of Clay Lacefield and Randy Bruno, Columbia University. The download_demo command will automatically download the file and store it in your caiman_data folder the first time you run it. To use the demo in your own dataset you can set:
fnames = [/path/to/file(s)].
End of explanation
display_movie = False
if display_movie:
    m_orig = cm.load_movie_chain(fnames)
    ds_ratio = 0.25
    m_orig.resize(1, 1, ds_ratio).play(
        q_max=99.5, fr=30, magnification=4)
Explanation: Viewing the data
Prior to analysis you can view the data by setting display_movie = False. The example file provided here is small, with FOV size 128 x 128. To view it larger we set magnification = 4 in the play command. This will result to a movie of resolution 512 x 512. For larger files, you will want to reduce this magnification factor (for memory and display purposes).
End of explanation
# dataset dependent parameters
fr = 30                             # imaging rate in frames per second
decay_time = 0.5
# motion correction parameters
strides = (48, 48)          # start a new patch for pw-rigid motion correction every x pixels
overlaps = (24, 24)         # overlap between pathes (size of patch strides+overlaps)
max_shifts = (6, 6)         # maximum allowed rigid shifts (in pixels)
max_deviation_rigid = 3     # maximum shifts deviation allowed for patch with respect to rigid shifts
pw_rigid = True             # flag for performing non-rigid motion correction
# parameters for source extraction and deconvolution
p = 0                       # order of the autoregressive system
gnb = 2                     # number of global background components
merge_thr = 0.75            # merging threshold, max correlation allowed
rf =  None                  # half-size of the patches in pixels. e.g., if rf=25, patches are 50x50
stride_cnmf = None          # amount of overlap between the patches in pixels
K = 60                      # number of components per patch
method_init = 'graph_nmf'   # initialization method (if analyzing dendritic data using 'sparse_nmf')
ssub = 1                    # spatial subsampling during initialization
tsub = 1                    # temporal subsampling during intialization
opts_dict = {'fnames': fnames,
            'fr': fr,
            'decay_time': decay_time,
            'strides': strides,
            'overlaps': overlaps,
            'max_shifts': max_shifts,
            'max_deviation_rigid': max_deviation_rigid,
            'pw_rigid': pw_rigid,
            'p': p,
            'nb': gnb,
            'rf': rf,
            'K': K, 
            'stride': stride_cnmf,
            'method_init': method_init,
            'rolling_sum': True,
            'only_init': True,
            'ssub': ssub,
            'tsub': tsub,
            'merge_thr': merge_thr}
opts = params.CNMFParams(params_dict=opts_dict)
#%% 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)
Explanation: Parameter setting
First we'll perform motion correction followed by CNMF.
One of the main differences when analyzing dendritic data, is that dendritic segments are not localized and can traverse significant parts of the FOV. They can also be spatially non-contiguous. To capture this property, CaImAn has two methods for initializing CNMF: sparse_nmf where a sparse non-negative matrix factorization is deployed, and graph-nmf, where the Laplacian of the pixel affinity matrix is used as a regularizer to promote spatial components that capture pixels that have similar timecourses. These can be selected with the parameter method_init.
In this demo we use the graph_nmf initialization method although sparse_nmf can also be used.
Since the dataset is small we can run CNMF without splitting the FOV in patches. This can be helpful because of the non-localized nature of the dendritic segments. To do that we set rf = None. Using patches is also supported as demonstrated below. Also note, that spatial and temporal downsampling can also be used to reduce the data dimensionality by appropriately modifying the parameters ssub and tsub. Here they are set to 1 since it is not necessary.
The graph NMF method was originally presented in the paper:
Cai, D., He, X., Han, J., & Huang, T. S. (2010). Graph regularized nonnegative matrix factorization for data representation. IEEE transactions on pattern analysis and machine intelligence, 33(8), 1548-1560.
End of explanation
mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion'))
%%capture
#%% Run piecewise-rigid motion correction using NoRMCorre
mc.motion_correct(save_movie=True)
m_els = cm.load(mc.fname_tot_els)
border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0 
    # maximum shift to be used for trimming against NaNs
Explanation: Motion correction
First we create a motion correction object and then perform the registration. The provided dataset has already been registered with a different method so we do not expect to see a lot of motion.
End of explanation
#%% compare with original movie
display_movie = False
if display_movie:
    m_orig = cm.load_movie_chain(fnames)
    ds_ratio = 0.2
    cm.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie,
                    m_els.resize(1, 1, ds_ratio)], 
                   axis=2).play(fr=60, q_max=99.5, magnification=4, offset=0)  # press q to exit
Explanation: Display registered movie next to the original
End of explanation
#%% MEMORY MAPPING
# memory map the file in order 'C'
fname_new = cm.save_memmap(mc.mmap_file, base_name='memmap_', order='C',
                           border_to_0=border_to_0) # exclude borders
# now load the file
Yr, dims, T = cm.load_memmap(fname_new)
images = np.reshape(Yr.T, [T] + list(dims), order='F') 
    #load frames in python format (T x X x Y)
#%% restart cluster to clean up memory
cm.stop_server(dview=dview)
c, dview, n_processes = cm.cluster.setup_cluster(
    backend='local', n_processes=None, single_thread=False)
Explanation: Memory mapping
End of explanation
cnm = cnmf.CNMF(n_processes, params=opts, dview=dview)
cnm = cnm.fit(images)
Explanation: Now run CNMF with the chosen initialization method
Note that in the parameter setting we have set only_init = True so only the initialization will be run.
End of explanation
CI = cm.local_correlations(images[::1].transpose(1,2,0))
CI[np.isnan(CI)] = 0
cnm.estimates.hv_view_components(img=CI)
Explanation: Display components
And compute the correlation image as a reference image
End of explanation
A = cnm.estimates.A.toarray().reshape(opts.data['dims'] + (-1,), order='F').transpose([2, 0, 1])
Nc = A.shape[0]
grid_shape = (np.ceil(np.sqrt(Nc/2)).astype(int), np.ceil(np.sqrt(Nc*2)).astype(int))
plt.figure(figsize=np.array(grid_shape[::-1])*1.5)
plt.imshow(montage(A, rescale_intensity=True, grid_shape=grid_shape))
plt.title('Montage of found spatial components');
plt.axis('off');
Explanation: Plot a montage array of all the components
We can also plot a montage of all the identified spatial footprints.
End of explanation
cnm2 = cnm.refit(images, dview=dview)
Explanation: Run CNMF to sparsify the components
The components found by the initialization have captured a lot of structure but they are not sparse. Using the refit command will sparsify them.
End of explanation
cnm2.estimates.hv_view_components(img=CI)
A2 = cnm2.estimates.A.toarray().reshape(opts.data['dims'] + (-1,), order='F').transpose([2, 0, 1])
Nc = A2.shape[0]
grid_shape = (np.ceil(np.sqrt(Nc/2)).astype(int), np.ceil(np.sqrt(Nc*2)).astype(int))
plt.figure(figsize=np.array(grid_shape[::-1])*1.5)
plt.imshow(montage(A2, rescale_intensity=True, grid_shape=grid_shape))
plt.title('Montage of found spatial components');
plt.axis('off');
Explanation: Now plot the components again
End of explanation
cnm2.estimates.play_movie(images, q_max=99.9,
                          magnification=4,
                          include_bck=False)
# press q to exit
Explanation: The resulting components are significantly more sparse and capture a lot of the data structure.
Play a movie with the results
The movie will show three panels with the original, denoised and residual movies respectively. 
The flag include_bck = False will remove the estimated background from the original and denoised panels.
To include it set it to True.
End of explanation
#cnm2.estimates.Cn = CI # save the correlation image for displaying in the background
#cnm2.save('dendritic_analysis.hdf5')
Explanation: Save the object
You can save the results of the analysis for future use and/or to load them on the GUI.
End of explanation
#%% restart cluster to clean up memory
cm.stop_server(dview=dview)
c, dview, n_processes = cm.cluster.setup_cluster(
    backend='local', n_processes=None, single_thread=False)
Explanation: Now try CNMF with patches
As mentioned above, the file shown here is small enough to be processed at once. For larger files patches and/or downsampling can be used. We show how to use patches below. Note that unlike in somatic imaging, for dendritic imaging we want the patches to be large and have significant overlap so enough structure of the dendritic segments is captured in the overlapping region. This will help with merging the components. However, a too large patch size can result to memory issues due to large amount of data loaded concurrently onto memory.
End of explanation
rf = 48   # size of each patch is (2*rf, 2*rf) 
stride_cnmf = 16
Kp = 25   # reduce the number of component as it is now per patch
opts.change_params({'rf': rf,
                    'stride': stride_cnmf,
                    'K': Kp});
cnm_patches = cnmf.CNMF(n_processes, params=opts, dview=dview)
cnm_patches = cnm_patches.fit(images)
cnm_patches.estimates.hv_view_components(img=CI)
Ap = cnm_patches.estimates.A.toarray().reshape(opts.data['dims'] + (-1,), order='F').transpose([2, 0, 1])
Nc = Ap.shape[0]
grid_shape = (np.ceil(np.sqrt(Nc/2)).astype(int), np.ceil(np.sqrt(Nc*2)).astype(int))
plt.figure(figsize=np.array(grid_shape[::-1])*1.5)
plt.imshow(montage(Ap, rescale_intensity=True, grid_shape=grid_shape))
plt.title('Montage of found spatial components');
plt.axis('off');
Explanation: Change some parameters to enable patches
To enable patches the rf parameter should be changed from None to some integer value. Moreover the amount of overlap controlled by stride_cnmf should also change. Finally, since the parameter K specifies components per patch it should be reduced compared to its value when operating on the whole FOV at once. It is important to experiment with different values for these parameters since results can sensitive, due to the non stereotypical shapes and activity patterns of dendrites.
End of explanation
cnm_patches2 = cnm_patches.refit(images, dview=dview)
Ap2 = cnm_patches2.estimates.A.toarray().reshape(opts.data['dims'] + (-1,), order='F').transpose([2, 0, 1])
Nc = Ap2.shape[0]
grid_shape = (np.ceil(np.sqrt(Nc/2)).astype(int), np.ceil(np.sqrt(Nc*2)).astype(int))
plt.figure(figsize=np.array(grid_shape[::-1])*1.5)
plt.imshow(montage(Ap2, rescale_intensity=True, grid_shape=grid_shape))
plt.title('Montage of found spatial components');
plt.axis('off');
cnm_patches2.estimates.play_movie(images, q_max=99.9,
                                  magnification=4,
                                  include_bck=False)
# press q to exit
Explanation: Refit to sparsify
As before we can pass the results through the CNMF refit function that can yield more sparse components.
End of explanation
AA = np.corrcoef(cnm2.estimates.A.toarray(), cnm_patches2.estimates.A.toarray(), rowvar=False)
plt.imshow(AA[:cnm2.estimates.A.shape[1], cnm2.estimates.A.shape[1]:])
plt.colorbar();
plt.xlabel('with patches')
plt.ylabel('without patches')
plt.title('Correlation coefficients');
Explanation: Compare the two approaches
To identify how similar components and the two approaches yield, we can compute the cross-correlation matrix between the spatial footprints derived from the two approaches.
End of explanation
#%% STOP CLUSTER and clean up log files
cm.stop_server(dview=dview)
log_files = glob.glob('*_LOG_*')
for log_file in log_files:
    os.remove(log_file)
Explanation: Stop the cluster when you're done
End of explanation
<END_TASK>