shape issue in update_spatial_components

[mschachter]

• Tell us about your operating system (Linux/macOS/Windows), Python version (2.7/3.x), and working environment (Spyder/Jupyter Notebook/other) OS: CENTOS 7 Python: 3.6 (anaconda 3) Environment: caiman environment from environment.yml Dataset: download here (67MB): Branch: dev

https://drive.google.com/file/d/14RVBCdqf-DSUvvRX-3tGYOHjSiotyHkT/view?usp=sharing

• Describe the issue that you are experiencing The movie (see link above) is a patch from a larger movie, and I don't think the patch has any cells. When I run caiman in one-photon mode using the code below, there is a failure that involves shapes. I spent a good amount of time investigating, and found some weird code in update_spatial_components that looks like this (line 234 of spatial.py):

    ff = np.where(np.sum(A_, axis=0) == 0)  # remove empty components
    if np.size(ff) > 0:
        ff = ff[0]
        print('eliminating {} empty spatial components'.format(len(ff)))
        A_ = np.delete(A_, list(ff[ff < nr]), 1)
        C = np.delete(C, list(ff[ff < nr]), 0)
        nr = nr - len(ff[ff < nr])
        if low_rank_background:
            background_ff = list(filter(lambda i: i >= nb, ff - nr))
            f = np.delete(f, background_ff, 0)
        else:
            background_ff = list(filter(lambda i: i >= 0, ff - nr))
            f = np.delete(f, background_ff, 0)
            b_in = np.delete(b_in, background_ff, 1)

This code doesn't make much sense to me - the variable "nr" is supposed to represent the number of cells, but there are several issues:

1. For the variable A_, the first index is the number of pixels, not the number of cells. Yet the variable "ff" is used as if it represents the number of cells.
2. The number of cells "nr" seems conflated with the number of background components "nb" later in the code.
3. The argument to the function update_spatial_components named "low_rank_background" is always true by default, because the value of "low_rank_background" from CNMFParams is never passed in from cnmf.py.

• Code to reproduce error Use this code to reproduce the error, after downloading the dataset:

  
  import os

from caiman import load_memmap from caiman.source_extraction.cnmf import CNMF from caiman.components_evaluation import estimate_components_quality_auto from caiman.source_extraction.cnmf.params import CNMFParams

set the params

frame_rate = 10.0 decay_time = 0.400 gSig = (9, 9) min_SNR = 3

cnmfe_params = CNMFParams(k=25, gSiz=(13, 13), decay_time=decay_time, min_pnr=2, min_SNR=min_SNR, p=1, merge_thresh=0.8, ssub=1, tsub=2, p_ssub=1, p_tsub=1, method_init='corr_pnr', min_corr=0.80, ssub_B=2, ring_size_factor=1.5,
low_rank_background=None, nb_patch=-1, gnb=-1)

load the file

root_dir = '/path/to/dataset' mov_file = os.path.join(root_dir, 'recording_20160517_172653-shapebug1_d1_64_d2_64_d3_1_order_C_frames4298.mmap')

Yr, dims, T = load_memmap(mov_file) Y = Yr.T.reshape((T,) + dims, order='F') print('*** Y.shape=',Y.shape)

fit the data

num_procs = 1 cnmfe = CNMF(num_procs, dview=None, params=cnmfe_params) cnmfe.fit(Y)

- **This is the output:**

WARNING:root:gnb=-1, hence setting keys nb_patch and low_rank_background in group patch automatically. * Y.shape= (4298, 64, 64) WARNING:root:Changing key init_batch in group online from 200 to 4298 WARNING:root:Changing key medw in group spatial from None to (3, 3) WARNING:root:Changing key se in group spatial from None to [[1 1 1] [1 1 1] [1 1 1]] WARNING:root:Changing key ss in group spatial from None to [[1 1 1] [1 1 1] [1 1 1]] WARNING:root:Changing key n_pixels_per_process in group preprocess from None to 4096 WARNING:root:Changing key n_pixels_per_process in group spatial from None to 4096 checking if missing data Noise Normalization Spatial Downsampling 1-photon Roi Extraction... One photon initialization.. 0 neurons have been initialized In total, 25 neurons were initialized. Compute Background Update Spatial Initializing update of Spatial Components computing the distance indicators /home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/scipy/sparse/compressed.py:774: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient. SparseEfficiencyWarning) memmaping Updating Spatial Components using lasso lars thresholding components eliminating 25 empty spatial components Computing residuals --- 7.87296199798584 seconds --- Removing tempfiles created Update Temporal Generating residuals /home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/scipy/sparse/dia.py:300: RuntimeWarning: divide by zero encountered in remainder c = np.arange(num_rows, dtype=np.intc) - (offsets % max_dim)[:, None] entering the deconvolution stopping: overall temporal component not changing significantly Initialization again 0 neurons have been initialized In total, 25 neurons were initialized. Merge Components No neurons merged! Update Spatial Initializing update of Spatial Components computing the distance indicators memmaping Updating Spatial Components using lasso lars thresholding components eliminating 25 empty spatial components Computing residuals --- 7.616950988769531 seconds --- Removing tempfiles created Update Temporal Generating residuals entering the deconvolution stopping: overall temporal component not changing significantly Compute Background Again Merge Components No neurons merged! Update Spatial Initializing update of Spatial Components computing the distance indicators memmaping Updating Spatial Components using lasso lars thresholding components Computing residuals --- 1.639448881149292 seconds --- Removing tempfiles created Update Temporal Generating residuals entering the deconvolution stopping: overall temporal component not changing significantly Return full Background False WARNING:root:NOT setting value of key use_init in group spatial, because no prior key existed... Initializing update of Spatial Components computing the distance indicators memmaping Updating Spatial Components using lasso lars thresholding components eliminating 4096 empty spatial components Computing residuals Traceback (most recent call last): File "scripts/zero_cells_bug.py", line 53, in cnmfe.fit(Y) File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/cnmf.py", line 480, in fit self.update_spatial(Yr, use_init=True) File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/cnmf.py", line 881, in update_spatial sn=self.estimates.sn, dims=self.dims, self.params.get_group('spatial')) File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/spatial.py", line 267, in update_spatial_components b = HALS4shape_bckgrnd(Y_resf, b_in, f, ind_b) File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/spatial.py", line 299, in HALS4shape_bckgrnd ((U[m, ind_pixels] - V[m].dot(B[ind_pixels].T)) / IndexError: index 0 is out of bounds for axis 0 with size 0

[agiovann]

@j-friedrich I had a look and the problem seems to be that the update_spatial function receives empty b and f components from the one-photon initialization, which then leads to a failure. I have set up the example in my desktop if you want to check it. I am not sure whether I should fix the function update_spatial, if this is a problem of parameters, or if it should be fixed in the initialization of 1-photon. Let me know what you think and I will proceed.

[j-friedrich]

Merely due to a wrong settings of parameters. cnmf.params.patch['only_init'] has to be True for 1p data.

[pgunn]

@j-friedrich It probably makes sense to have the function check for this and inform the user.

[j-friedrich]

Yes, could maybe also set it automatically to True if center_psf is True, ring_size_factor is not None and method_init=='corr_pnr', which together specify 1p data.

[mschachter]

Hi @j-friedrich, setting only_patch_init fixes the immediate error, but there is a new error that pops up, based on the shape of either A or C.

Here is the code to reproduce the problem:

import os

from caiman import load_memmap
from caiman.source_extraction.cnmf import CNMF
from caiman.components_evaluation import estimate_components_quality_auto
from caiman.source_extraction.cnmf.params import CNMFParams

# set the params
frame_rate = 10.0
decay_time = 0.400
gSig = (9, 9)
min_SNR = 3

cnmfe_params = CNMFParams(k=25, gSiz=(13, 13), decay_time=decay_time, min_pnr=2, 
                          min_SNR=min_SNR, p=1, merge_thresh=0.8,
                          ssub=1, tsub=2, p_ssub=1, p_tsub=1, method_init='corr_pnr',
                          min_corr=0.80, ssub_B=2, ring_size_factor=1.5,                          
                          low_rank_background=None, nb_patch=-1, gnb=-1,
                          only_init_patch=True)

# load the file
root_dir = '/path/to/dataset'
mov_file = os.path.join(root_dir, 'recording_20160517_172653-shapebug1_d1_64_d2_64_d3_1_order_C_frames_4298_.mmap')

Yr, dims, T = load_memmap(mov_file)
Y = Yr.T.reshape((T,) + dims, order='F')
print('*** Y.shape=',Y.shape)

# fit the data
num_procs = 1
cnmfe = CNMF(num_procs, dview=None, params=cnmfe_params)
cnmfe.fit(Y)

Here is the new error:

Traceback (most recent call last):
  File "scripts/zero_cells_bug.py", line 53, in <module>
    cnmfe.fit(Y)
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/cnmf.py", line 448, in fit
    self.compute_residuals(Yr)
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/cnmf.py", line 663, in compute_residuals
    self.estimates.YrA = (YA - (AA.T.dot(Cf)).T)[:, :self.estimates.A.shape[-1]].T
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/scipy/sparse/base.py", line 302, in dot
    return self * other
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/scipy/sparse/base.py", line 405, in __mul__
    raise ValueError('dimension mismatch')
ValueError: dimension mismatch

[j-friedrich]

Another parameter. For 1p data cnmf.params.init['center_psf'] is supposed to be True. We actually treat center_psf almost like a flag for whether the data is 1p or 2p, though one could use centering also for 2p and set ring_size_factor to None. There's however indeed and error in line 447 of cnmf.py. Should read if not (self.params.get('init', 'method_init') == 'corr_pnr' and self.params.get('init', 'ring_size_factor') is not None): Your script could be an exotic but reasonable choice for using the ring model for 2p data, although one would usually use a low rank estimate. 2p usually doesn't require center_psf=True in order to roughly remove the background when finding and initializing neurons. With the above change in line 447 your script doesn't run through yet, will sort this out later.

Ropening until we added some parameter checks and inform the user.

[mschachter]

Thanks @j-friedrich . I can confirm that there is a new error when I set center_psf to True.

It was not my intention to write such an exotic script :) I was trying to follow your instructions to use the ring model with full background on 1P data, but definitely I realize that center_psf should have been set to True now.

Traceback (most recent call last):
  File "scripts/zero_cells_bug.py", line 53, in <module>
    cnmfe.fit(Y)
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/cnmf.py", line 475, in fit
    self.estimates.normalize_components()
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/caiman/source_extraction/cnmf/estimates.py", line 598, in normalize_components
    self.f = nB_mat * self.f
  File "/home/mschachter/anaconda3/envs/caiman/lib/python3.6/site-packages/scipy/sparse/base.py", line 405, in __mul__
    raise ValueError('dimension mismatch')
ValueError: dimension mismatch

[j-friedrich]

Fixed the issue in branch issue362, so that both of your scripts run through, even the one with only_init_patch=False. The scenario would be to folllow up the standard pipeline, that's run within the initialization using the ring model, by few spatial and temporal updates while keeping the full-rank background fixed. Remains to see whether that slightly improves results.

[mschachter]

Thanks again @j-friedrich! I can confirm that the issue362 branch fixes my issue with that patch, and does not seem to do anything weird to previously working patches. I can now continue on with my benchmarking of 1P CNMFe, I'll keep an eye out for when issue362 is merged with dev.

[epnev]

@mschachter Is this related to the discussion you had with @j-friedrich in slack? (Just returned from vacation and I'm trying to catch up)

[mschachter]

He definitely made changes to the code that works with the variable background_ff, so I think it is related.

[epnev]

@mschachter I believe this has been fixed on dev with 709fe653c65298ff81549a4e5268a8979d8b5b9c so close if that's the case.

[j-friedrich]

Yes, fixed in dev. Not in master, but fine with closing.