Merge pull request #863 from flatironinstitute/dev

Merge dev to master for release

caiman/base/rois.py

@@ -191,7 +191,8 @@ def nf_match_neurons_in_binary_masks(masks_gt,
                                      labels=['Session 1', 'Session 2'],
                                      cmap='gray',
                                      D=None,
-                                     enclosed_thr=None):
+                                     enclosed_thr=None,
+                                     colors=['red', 'white']):
     """
     Match neurons expressed as binary masks. Uses Hungarian matching algorithm
 
@@ -296,27 +297,27 @@ def nf_match_neurons_in_binary_masks(masks_gt,
             font = {'family': 'Myriad Pro', 'weight': 'regular', 'size': 10}
             pl.rc('font', **font)
             lp, hp = np.nanpercentile(Cn, [5, 95])
-            ses_1 = mpatches.Patch(color='red', label=labels[0])
-            ses_2 = mpatches.Patch(color='white', label=labels[1])
+            ses_1 = mpatches.Patch(color=colors[0], label=labels[0])
+            ses_2 = mpatches.Patch(color=colors[1], label=labels[1])
             pl.subplot(1, 2, 1)
             pl.imshow(Cn, vmin=lp, vmax=hp, cmap=cmap)
-            [pl.contour(norm_nrg(mm), levels=[level], colors='w', linewidths=1) for mm in masks_comp[idx_tp_comp]]
-            [pl.contour(norm_nrg(mm), levels=[level], colors='r', linewidths=1) for mm in masks_gt[idx_tp_gt]]
+            [pl.contour(norm_nrg(mm), levels=[level], colors=colors[1], linewidths=1) for mm in masks_comp[idx_tp_comp]]
+            [pl.contour(norm_nrg(mm), levels=[level], colors=colors[0], linewidths=1) for mm in masks_gt[idx_tp_gt]]
             if labels is None:
                 pl.title('MATCHES')
             else:
-                pl.title('MATCHES: ' + labels[1] + '(w), ' + labels[0] + '(r)')
+                pl.title('MATCHES: ' + labels[1] + f'({colors[1][0]}), ' + labels[0] + f'({colors[0][0]})')
             pl.legend(handles=[ses_1, ses_2])
             pl.show()
             pl.axis('off')
             pl.subplot(1, 2, 2)
             pl.imshow(Cn, vmin=lp, vmax=hp, cmap=cmap)
-            [pl.contour(norm_nrg(mm), levels=[level], colors='w', linewidths=1) for mm in masks_comp[idx_fp_comp]]
-            [pl.contour(norm_nrg(mm), levels=[level], colors='r', linewidths=1) for mm in masks_gt[idx_fn_gt]]
+            [pl.contour(norm_nrg(mm), levels=[level], colors=colors[1], linewidths=1) for mm in masks_comp[idx_fp_comp]]
+            [pl.contour(norm_nrg(mm), levels=[level], colors=colors[0], linewidths=1) for mm in masks_gt[idx_fn_gt]]
             if labels is None:
-                pl.title('FALSE POSITIVE (w), FALSE NEGATIVE (r)')
+                pl.title(f'FALSE POSITIVE ({colors[1][0]}), FALSE NEGATIVE ({colors[0][0]})')
             else:
-                pl.title(labels[1] + '(w), ' + labels[0] + '(r)')
+                pl.title(labels[1] + f'({colors[1][0]}), ' + labels[0] + f'({colors[0][0]})')
             pl.legend(handles=[ses_1, ses_2])
             pl.show()
             pl.axis('off')

caiman/motion_correction.py

@@ -179,12 +179,15 @@ def __init__(self, fname, min_mov=None, dview=None, max_shifts=(6, 6), niter_rig
         """
         if 'ndarray' in str(type(fname)):
             logging.info('Creating file for motion correction "tmp_mov_mot_corr.hdf5"')
-            cm.movie(fname).save('./tmp_mov_mot_corr.hdf5')
+            cm.movie(fname).save('./tmp_mov_mot_corr.hdf5') # FIXME don't write to the current directory!
             fname = ['./tmp_mov_mot_corr.hdf5']
 
         if type(fname) is not list:
             fname = [fname]
 
+        if type(gSig_filt) is tuple:
+            gSig_filt = list(gSig_filt) # There are some serializers down the line that choke otherwise
+
         self.fname = fname
         self.dview = dview
         self.max_shifts = max_shifts

caiman/source_extraction/cnmf/utilities.py

@@ -968,10 +968,10 @@ def get_file_size(file_name, var_name_hdf5='mov'):
                 if loading from hdf5 name of the dataset to load
 
         Returns:
-            dims: list
+            dims: tuple
                 dimensions of FOV
 
-            T: list
+            T: int or tuple of int
                 number of timesteps in each file
     """
     if isinstance(file_name, pathlib.Path):
@@ -1059,6 +1059,10 @@ def get_file_size(file_name, var_name_hdf5='mov'):
         else:
             dims, T = zip(*[get_file_size(fn, var_name_hdf5=var_name_hdf5)
                 for fn in file_name])
+            if len(set(dims)) > 1:
+                raise Exception('Files have FOVs with different sizes')
+            else:
+                dims = dims[0]
     else:
         raise Exception('Unknown input type')
     return dims, T

caiman/source_extraction/volpy/spikepursuit.py

@@ -145,12 +145,9 @@ def volspike(pars):
                     
                 thresh: float
                     threshold of the signal
-
-                spatial_filter: 2-d array
-                    spatial filter of the neuron in the FOV
                     
                 weights: 2-d array
-                    ridge regression coefficients for fitting reconstructed signal
+                    ridge regression coefficients for fitting reconstructed signal 
                 
                 locality: boolean
                     False if the maximum of spatial filter is not in the initial ROI
@@ -226,12 +223,17 @@ def volspike(pars):
     # compute the initial trace
     if weights_init is None:
         t0 = np.nanmean(data_hp[:, bw.ravel()], 1)
-    else:
-        t0 = np.matmul(data_hp, weights_init[1:])  
+    else: 
+        print('Reuse weights')
+        t0 = np.matmul(data_hp, weights_init.ravel()) # reuse weights
     t0 = t0 - np.mean(t0)
 
     # remove any variance in trace that can be predicted from the background principal components
-    Ub, Sb, Vb = svds(data_hp[:, notbw.ravel()], args['nPC_bg'])
+    data_svd = data_hp[:, notbw.ravel()]
+    if data_svd.shape[1] < args['nPC_bg'] + 1:
+        raise Exception(f'Too few pixels ({data_svd.shape[1]}) for background extraction (at least {args["nPC_bg"]} needed);'
+                        f'please decrease context_size and censor_size')
+    Ub, Sb, Vb = svds(data_svd, args['nPC_bg'])
     alpha = args['nPC_bg'] * args['ridge_bg']    # square of F-norm of Ub is equal to number of principal components
     reg = Ridge(alpha=alpha, fit_intercept=False, solver='lsqr').fit(Ub, t0)
     t0 = np.double(t0 - np.matmul(Ub, reg.coef_))
@@ -246,7 +248,10 @@ def volspike(pars):
     output['rawROI']['t'] = t0.copy()
     output['rawROI']['ts'] = ts.copy()
     output['rawROI']['spikes'] = spikes.copy()
-    output['rawROI']['spatial_filter'] = bw.copy()
+    if weights_init is None:
+        output['rawROI']['weights'] = bw.copy()
+    else:
+        output['rawROI']['weights'] = weights_init.copy()
     output['rawROI']['t'] = output['rawROI']['t'] * np.mean(t0[output['rawROI']['spikes']]) / np.mean(
         output['rawROI']['t'][output['rawROI']['spikes']])  # correct shrinkage
     output['rawROI']['templates'] = templates
@@ -284,14 +289,14 @@ def volspike(pars):
                             kernel_std_x=sigma, kernel_std_y=sigma,
                             borderType=cv2.BORDER_REPLICATE), data_hp.shape)))
 
-    # refine spatial filter and estimate spike times for several iterations 
+    # refine weights and estimate spike times for several iterations 
     for iteration in range(args['n_iter']):
         if iteration == args['n_iter'] - 1:
             do_plot = args['do_plot']
         else:
             do_plot = False
 
-        # update spatial filter
+        # update weights
         tr = np.single(t_rec.copy())
         if args['weight_update'] == 'NMF':
             C = np.array([tr, np.ones_like(tr)])  # constant baselines as 2nd component
@@ -310,15 +315,6 @@ def volspike(pars):
             weights = Ri.coef_
             weights[0] = Ri.intercept_
 
-        spatial_filter = np.empty_like(weights)
-        spatial_filter[:] = weights
-        spatial_filter = movie.gaussian_blur_2D(np.reshape(spatial_filter[1:],
-                                                          ref.shape, order='C')[np.newaxis, :, :],
-                                               kernel_size_x=np.int(2 * np.ceil(2 * sigma) + 1),
-                                               kernel_size_y=np.int(2 * np.ceil(2 * sigma) + 1),
-                                               kernel_std_x=sigma, kernel_std_y=sigma,
-                                               borderType=cv2.BORDER_REPLICATE)[0]
-
         # update the signal            
         t = np.matmul(recon, weights)
         t = t - np.mean(t)
@@ -328,6 +324,7 @@ def volspike(pars):
         t = t - np.matmul(Ub, b)
 
         # correct shrinkage
+        weights = weights * np.mean(t0[spikes]) / np.mean(t[spikes])
         t = np.double(t * np.mean(t0[spikes]) / np.mean(t[spikes]))
 
         # estimate spike times
@@ -362,14 +359,11 @@ def volspike(pars):
     else:
         locality = True
 
-    # spatial filter in the FOV
+    # weights in the FOV
     weights = np.reshape(weights[1:],ref.shape, order='C')
     weights_FOV = np.zeros(images.shape[1:])
     weights_FOV[Xinds[0]:Xinds[-1] + 1, Yinds[0]:Yinds[-1] + 1] = weights
 
-    spatial = np.zeros(images.shape[1:])
-    spatial[Xinds[0]:Xinds[-1] + 1, Yinds[0]:Yinds[-1] + 1] = spatial_filter
-
     # subthreshold activity extraction    
     t_sub = t.copy() - t_rec
     t_sub = signal_filter(t_sub, args['sub_freq'], fr, order=5, mode='low') 
@@ -386,7 +380,6 @@ def volspike(pars):
     output['templates'] = templates
     output['snr'] = snr
     output['thresh'] = thresh
-    output['spatial_filter'] = spatial    
     output['weights'] = weights_FOV
     output['locality'] = locality    
     output['context_coord'] = np.transpose(np.vstack((Xinds[[0, -1]], Yinds[[0, -1]])))

caiman/source_extraction/volpy/utils.py

@@ -274,4 +274,5 @@ def update(val):
     s_comp.on_changed(update)
     s_comp.set_val(0)
     fig.canvas.mpl_connect('key_release_event', arrow_key_image_control)
-    plt.show()
+    plt.show()
+

demos/general/demo_pipeline_voltage_imaging.py

@@ -183,8 +183,6 @@ def main():
     fig, axs = plt.subplots(1, 2)
     axs[0].imshow(summary_images[0]); axs[1].imshow(ROIs.sum(0))
     axs[0].set_title('mean image'); axs[1].set_title('masks')
-
-
 
 # %% restart cluster to clean up memory
     cm.stop_server(dview=dview)
@@ -194,7 +192,7 @@ def main():
 # %% parameters for trace denoising and spike extraction
     ROIs = ROIs                                   # region of interests
     index = list(range(len(ROIs)))                # index of neurons
-    weights = None                                # reuse spatial weights 
+    weights = None                                # if None, use ROIs for initialization; to reuse weights check reuse weights block 
 
     context_size = 35                             # number of pixels surrounding the ROI to censor from the background PCA
     visualize_ROI = False                         # whether to visualize the region of interest inside the context region
@@ -235,7 +233,7 @@ def main():
 #%% TRACE DENOISING AND SPIKE DETECTION
     vpy = VOLPY(n_processes=n_processes, dview=dview, params=opts)
     vpy.fit(n_processes=n_processes, dview=dview)
-
+   
 #%% visualization
     display_images = True
     if display_images:
@@ -248,7 +246,7 @@ def main():
     if display_images:
         mv_all = utils.reconstructed_movie(vpy.estimates.copy(), fnames=mc.mmap_file,
                                            idx=idx, scope=(0,1000), flip_signal=flip_signal)
-        mv_all.play(fr=40)
+        mv_all.play(fr=40, magnification=3)
 
 #%% save the result in .npy format 
     save_result = True
@@ -257,6 +255,16 @@ def main():
         vpy.estimates['params'] = opts
         save_name = f'volpy_{os.path.split(fnames)[1][:-5]}_{threshold_method}'
         np.save(os.path.join(file_dir, save_name), vpy.estimates)
+
+#%% reuse weights 
+# set weights = reuse_weights in opts_dict dictionary
+    estimates = np.load(os.path.join(file_dir, save_name+'.npy'), allow_pickle=True).item()
+    reuse_weights = []
+    for idx in range(ROIs.shape[0]):
+        coord = estimates['context_coord'][idx]
+        w = estimates['weights'][idx][coord[0][0]:coord[1][0]+1, coord[0][1]:coord[1][1]+1] 
+        plt.figure(); plt.imshow(w);plt.colorbar(); plt.show()
+        reuse_weights.append(w)
 
 # %% STOP CLUSTER and clean up log files
     cm.stop_server(dview=dview)

demos/notebooks/demo_dendritic.ipynb

@@ -382,6 +382,81 @@
     "plt.axis('off');"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Component evaluation\n",
+    "the components are evaluated in three ways:\n",
+    "   a) the shape of each component must be correlated with the data\n",
+    "   b) a minimum peak SNR is required over the length of a transient\n",
+    "   c) each shape passes a CNN based classifier"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "min_SNR = 4  # signal to noise ratio for accepting a component\n",
+    "rval_thr = 0.85  # space correlation threshold for accepting a component\n",
+    "cnn_thr = 0.99  # threshold for CNN based classifier\n",
+    "cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected\n",
+    "\n",
+    "cnm2.params.set('quality', {'decay_time': decay_time,\n",
+    "                           'min_SNR': min_SNR,\n",
+    "                           'rval_thr': rval_thr,\n",
+    "                           'use_cnn': False, # do not use cnn for dentritic data\n",
+    "                           'min_cnn_thr': cnn_thr,\n",
+    "                           'cnn_lowest': cnn_lowest})\n",
+    "\n",
+    "cnm2.estimates.evaluate_components(images, cnm2.params, dview=dview)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot contours of selected and rejected components"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# accepted components\n",
+    "cnm2.estimates.hv_view_components(img=CI, idx=cnm2.estimates.idx_components)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# rejected components\n",
+    "cnm2.estimates.hv_view_components(img=CI, idx=cnm2.estimates.idx_components_bad)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Extract DF/F values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -426,8 +501,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#cnm2.estimates.Cn = CI # save the correlation image for displaying in the background\n",
-    "#cnm2.save('dendritic_analysis.hdf5')"
+    "cnm2.estimates.Cn = CI # save the correlation image for displaying in the background\n",
+    "cnm2.save('dendritic_analysis.hdf5')"
    ]
   },
   {
@@ -610,9 +685,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [conda env:caiman-nwb] *",
+   "display_name": "Python 3",
    "language": "python",
-   "name": "conda-env-caiman-nwb-py"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -624,7 +699,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,