zhoupc / cnmf_e Goto Github PK

Hello. I am getting the error below from neuron.update_temporal_parallel(use_parallel) when running demo_large_data_1p.m.

I went back and tried the unchanged script with the demo data and did not get this error, so it seems possible my parameters (or data) are causing this error downstream. Would you have any idea what may have caused this and how I may correct it?

Thanks for the help, it is much appreciated.

The error:

`Error using foopsi_oasisAR1>update_g (line 136)
Index in position 2 exceeds array bounds.

Error in foopsi_oasisAR1 (line 100)
[solution, active_set, g, spks] = update_g(y-b, active_set,lam, smin);

Error in deconvolveCa (line 119)
[c, s, options.b, options.pars] = foopsi_oasisAR1(y-options.b, options.pars, options.lambda, ...

Error in HALS_temporal (line 92)
[ck, sk, tmp_options]= deconvolveCa(ck_raw, deconv_options, 'maxIter', 20, 'sn', sn_psd, 'pars', kernel_pars{k});

Error in Sources2D/update_temporal_parallel (line 112)
parfor mpatch=1:(nr_patch*nc_patch)`

Hi, From Wiki and the paper, it suggests S(i) is the inferred number of spike events at i-th frame. I am thinking to use S(i) as a proxy to calculate firing rate. However, in CalmAn (Python package implementation of CNMF), it is emphasized that S cannot be directly used as number of spikes or spiking rate, and instead, it is just loosely proportional to spiking rate. It suggests people to use dF_F (F_dff) to deconvolve and further set some amplitude threshold. I also have seen papers taking CNMFe results and treat each S(i) as a potential spiking event and ignore the amplitude while calculating firing rate (ref: https://www.biorxiv.org/content/10.1101/292953v3.full.pdf). I'm thus curious to know which is the right way to do? Treat each S(i) as one spike event, or use the number of S(i) as spiking rate? Thanks a lot!

I get an error about incorrect matrix dimensions in the demo_endoscope.m line 249 - number of column of neuron.A is different from the number of rows in matrix neuron.C This is most probably due to an error in the updateTemporal_endoscope.m where the matrix A is not updated. Probably, a line

obj.A = A(:,~ind_del);

is missing at the end of the function.

Hi I'm using the demo for large_data_1p. After initializing the residual, the script does update_spatial_parallel and I got a warning:

Warning: Matrix is singular to working precision.

In com (line 24)
In determine_search_location (line 58)
In Sources2D/update_spatial_parallel (line 66)

Then an error:

Error using determine_search_location (line 66)
Input matrix contains NaN or Inf.

Error in Sources2D/update_spatial_parallel (line 66)
IND = sparse(logical(determine_search_location(obj.A, search_method, options)));

I'm having trouble understanding why this error occurs and what it means. Does it have to do with the step where neurons are picked from the residual? The seeds that were picked in that step hardly coincide with the seeds picked during initialization.

@zhoupc

There is a bug at https://github.com/zhoupc/CNMF_E/blob/master/ca_source_extraction/%40Sources2D/Sources2D.m#L1810.

if ~exist('zip_file_path', 'var') || ~isempty(zip_file_path)

should be changed to

if exist('zip_file_path', 'var') && ~isempty(zip_file_path)

to avoid errors in the default scripts when the method is called without a zip file specified (e.g. neuron.save_workspace();). See https://github.com/zhoupc/CNMF_E/blob/master/demos/demo_large_data_1p.m#L215.

@zhoupc

Can you provide an option to either put the temporary folder created for the matfile in a user defined sub-folder or delete the folder after a successful run on CNMF-E? This folder can get large and take up unnecessary space after CNMF-E has been run.

e.g. in the below section of Sources2D, instead of [] for dir_output input, this should be a user modifiable parameter.
https://github.com/zhoupc/CNMF_E/blob/master/ca_source_extraction/%40Sources2D/Sources2D.m#L239

Hello, I encountered the following error when running the demo_batch_2p.m and demo_batch_1p.m files:
/*************************
Error using save
No variable matched the pattern 'save_*'.

Error in Sources2D/save_workspace (line 1781)
evalin('caller', sprintf('save(''%s'', ''neuron'', ''save_'', ''show_'', ''use_parallel'', ''with_*'', ''-v7.3''); ', file_path));

Error in Sources2D/save_workspace_batch (line 1821)
obj.save_workspace();

Error in demo_batch_2p (line 128)
neuron.save_workspace_batch();
*****************************/

I just ran the program without changing anything.

P.S. : I cloned the git repository as suggested in a previous issue, and the problem still existed.

Hello,

I'm a graduate student in Prof. Minwhan Jung's lab, KAIST, South Korea.
I have an issue that appears when running demo_large_data_1p.m, in the first update of spatial components, after
[center_res, Cn_res, PNR_res] =neuron.initComponents_residual_parallel([], save_initialization, use_parallel, min_corr_res, min_pnr_res, seed_method_res);.

The error message leads me to line 66,
IND = sparse(logical(determine_search_location(obj.A, search_method, options)));
and again to line 66 of determine_search_location, which is a parfor.
I run manually this parfor and found that in line 69, [V,D] = eig(Vr{i}); was stopping me, as a column of obj.A was a zero sparse. I realized that some zero-pixel neurons were existing, but the fixing of this false positives (remove_false_positives) is placed after this spatial update in the demo, so I'm opening this issue.

I tried to overcome with starting with higher min_pnr and min_corr but did not work.

Thanks.
Young Ju.

Hello

Thanks for sharing your code. I've been trying it out with the test dataset you provide (run demos/demo_large_data_2p.m) and I ran into the error below. Any idea what it could be? My best guess is that I'm maybe missing some of the dependencies.

Thanks for your help

Tiago

You have ran 1 initialization(s). 
*  1:	 LOGS_13-Dec_18_31_16
	-------------------------- GUIDE --------------------------
	use  i to select the one you want to continue your analysis
	type -i to remove some previous results from your hard drive
	type anything if you want to start a new initialization
	--------------------------  END  --------------------------
* make your choice:    1
* make your choice:    i
Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers.
10 neurons have been detected
10 neurons have been detected
20 neurons have been detected
20 neurons have been detected
Patch ( 1,  2) is done.  2 X  2 patches in total. 
Patch ( 2,  1) is done.  2 X  2 patches in total. 
10 neurons have been detected
10 neurons have been detected
20 neurons have been detected
Patch ( 2,  2) is done.  2 X  2 patches in total. 
Patch ( 1,  1) is done.  2 X  2 patches in total. 

-----------------UPDATE BACKGROUND---------------------------
Patch ( 1,  2) is done.  2 X  2 patches in total. 
Patch ( 2,  1) is done.  2 X  2 patches in total. 
Patch ( 1,  1) is done.  2 X  2 patches in total. 
Patch ( 2,  2) is done.  2 X  2 patches in total. 
Finished updating background using svd model.

---------------PICK NEURONS FROM THE RESIDUAL----------------
Patch ( 1,  2) is done.  2 X  2 patches in total. 
Patch ( 2,  1) is done.  2 X  2 patches in total. 
Patch ( 1,  1) is done.  2 X  2 patches in total. 
Patch ( 2,  2) is done.  2 X  2 patches in total. 

-----------------UPDATE TEMPORAL---------------------------
Undefined function or variable 'A_prev'.

Error in Sources2D/update_temporal_parallel (line 112)
    parfor mpatch=1:(nr_patch*nc_patch)

Error in demo_large_data_2p (line 131)
neuron.update_temporal_parallel(use_parallel);

Error in run (line 96)
evalin('caller', [script ';']);

@zhoupc
Related to unresolved #31.

Please add Y(isnan(Y)) = 0; right after https://github.com/zhoupc/CNMF_E/blob/master/ca_source_extraction/endoscope/distribute_data.m#L183.

This fixes CNMF-E crashing out in https://github.com/zhoupc/CNMF_E/blob/master/ca_source_extraction/endoscope/connectivity_constraint.m.

IMPORTANT NOTE ABOUT MATLAB R2017a or later: a number of users are reporting issues with
CVX and R2017a; and indeed, we’re seeing them too. At the moment, we must assume that MATLAB
R2017a is *not* compatible with CVX. If you need to use CVX, you will need to stick with R2016b or
earlier. We do not have a timeline for a fix yet.

Is there anyway around this? An alternative de-noising dependency we could use?

@zhoupc

The current codebase throws an error at
https://github.com/zhoupc/CNMF_E/blob/master/ca_source_extraction/endoscope/connectivity_constraint.m#L13
if the movie contains NaNs (as might be the case after motion correction when the borders are cropped with NaNs). Re-ran a couple of movies after setting the NaNs to zero and it ran fine without errors.

Is there a way to have CNMF-E automatically account for this when running movies (e.g. have a inputMovie(isnan(inputMovie)) = 0; % or mean pixel value) so we don't have to alter or duplicate existing datasets? e.g. updating `Sources2D.load_data' method? https://github.com/zhoupc/CNMF_E/blob/master/ca_source_extraction/%40Sources2D/Sources2D.m#L1073

Dear @zhoupc

I have an issue which I can't resolve. Can you provide me some help?

Thank you in advance!

Error using zeros
Size inputs must be integers.

Error in get_patch_data (line 70)
Y = zeros(nr, nc, nframes, 'like', cast(0,dtype));

Error in Sources2D/estimate_noise (line 366)
Ypatch = get_patch_data(mat_data, [r0, r1, c0, c1], frame_range, false);

Error in Sources2D/getReady (line 251)
obj.P.sn = obj.estimate_noise();

Hello! I am having an issue using the script to save my results and workspace.
When I run

%% save results
results = neuron.obj2struct();
eval(sprintf('save %s%s%s_results.mat results', dir_nm, filesep, file_nm));

I receive the following messages:
Error using save
Variable 'Expansion' not found.

Error in demo_endoscope (line 236)
eval(sprintf('save %s%s%s_results.mat results', dir_nm, filesep, file_nm));

When I run %% save the workspace for future analysis
neuron.save_workspace();

I receive
Reference to non-existent field 'log_folder'.

Error in Sources2D/save_workspace (line 1413)
file_path = [obj.P.log_folder, strrep(get_date(), ' ', '_'), '.mat'];

Error in demo_endoscope (line 242)
neuron.save_workspace();

When I try to run display contours of the neurons, I get the error below.

Error: File: demo_endoscope.m Line: 202 Column: 40
The construct "Ysignal(...end...)" is ambiguous in this context, because "Ysignal" cannot be ascertained to be either the name of a variable or of a function. To make it a variable, assign to it; to allow it to become a function at execution time, replace "end" with a call to LENGTH, SIZE, or NUMEL on the desired array.

Firstly, thank you for sharing this project. I noticed that after the last update the program crashes here:

Do you maybe mean patch_sz = [ceil(d1/x) ceil(d2/x)]; ?

-----------------UPDATE BACKGROUND---------------------------
Error using fit_ring_model (line 39)
Sparse single array arithmetic operations are not supported.

Error in Sources2D/update_background_parallel (line 176)
parfor mpatch=1:(nr_patch*nc_patch)

Error in demo_large_data_1p (line 139)
neuron.update_background_parallel(use_parallel);

Hello,

I think you want to compare with "spline" (the "p" is missing).

edit: Found another one

ind_succwss should be ind_success .

Dear @zhoupc ,

Do we have any function in the toolbox that can reconstruct the neuron.C_raw if I input neuron.A and the original .avi Video of the dataset?
My question is because I am making some changes to the video before loading it to CNMF-E. I get the neuron object with neuron.A and neuron.C_raw variables, but this C_raw is biased by the changes I have made to the video. I wanted to use the cells' footprints of neuron.A and the original video to reconstruct the raw calcium traces. Only making matrix multiplication doesn't seem to be enough, since the code seems to have a baseline subtraction, right? Is there any easy way to do it?

Thank you very much in advance!

Best regards,
Daniel

After running initialize A and C, the neuron.S is filled, but when I run the 'merge neurons' cell and click through the potential merges, I get the following error and then subsequently neuron.S is empty.

If I accept all the merge suggestions, the error does not happen.

Hi PC -

I've been encountering a sporadic error (throws this error consistently for only some .avi files):

Error using Sources2D/initComponents_parallel (line 310)
Size arguments must be real integers.

Error in CNMFEpipeline1 (line 128)
[center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);

Any ideas on how to bypass this? Thanks!

Hi,

Thank you for sharing this very useful package. I greatly appreciate all of the work you've put into this for batch processing, and parallel processing. Here's one possibly way to make it better. If you agree, I'll try to make a pull request.

The memory usage for the A matrices (A_new, A_, etc.) could be improved. These are the matrices that hold the spatial mask for each neuron. My data is particularly large with > 10k neurons and a high resolution ~ 2.5 megapixels, so the A matrices (which are neurons x pixels) can quickly take up several hundred GB if they aren't sparse, but they are only tens of MB if sparse. I can see that you've tried to use sparse matrices where possible, but since there is a need to switch from 2D to 3D, and 3D sparse matrices aren't allowed, full() is often called, which is where I run into an out of memory error. I found a ND sparse class on the file exchange (link below). It seems like a fantastic implementation of sparse arrays. I was able to make A ndSparse throughout the code and it seems to work okay (removes the need to call full() and reshape works fine). Can you think of any reason that would break other things?

https://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays

Should I make a pull request for this? Or do you see this as not needed or possibly breaking other things? If so, that's okay too.

Thanks,
Tom

CNMF_E Version: 1.1.1

Movie file: (3.1GB) https://drive.google.com/file/d/1kNDgHwgfn2bCtEjH5g7RFeIwihGdEGx7/view?usp=sharing

Non-default Parameters:
gSiz: [10, 10]
K: 37
patch_dims: [25, 25]
gSig: [0, 0]
min_pnr: 5
min_corr: 0.8
memory_size_to_use: 64
memory_size_per_patch: 8
max_tau: 0.400

Description: when I run CNMF_E on the movie listed above with the given parameters, it eventually fails during a deconvolution step with this error:

Error using foopsi_oasisAR1>update_g (line 144)
Index exceeds matrix dimensions.Error in foopsi_oasisAR1 (line 106)
           [solution, active_set, g, spks] = update_g(y-b, active_set,lam, smin, g_range);Error in deconvolveCa (line 118)
           [c, s, options.b, options.pars] = foopsi_oasisAR1(y-options.b, options.pars, options.lambda, ...Error in HALS_temporal (line 94)
           [ck, sk, tmp_options]= deconvolveCa(ck_raw, deconv_options, 'maxIter', 20, 'sn', sn_psd, 'pars', kernel_pars{k});Error in Sources2D/update_temporal_parallel (line 112)
   parfor mpatch=1:(nr_patch*nc_patch)Error in isx.cnmfe.run_cnmfe (line 261)
   neuron.update_temporal_parallel(use_parallel);

I debugged this issue, and found that size(active_set, 1) == 0 for one of the traces. I fixed this problem by modifying the function of foopsi_oasisAR1 in one place (look for the string %%%%%%%%%% to see changes):

function [c, s, b, g, active_set] = foopsi_oasisAR1(y, g, lam, smin, optimize_b,...
    optimize_g, decimate, maxIter, tau_range, gmax)
%% Infer the most likely discretized spike train underlying an AR(1) fluorescence trace
% Solves the sparse non-negative deconvolution problem
%  min 1/2|c-y|^2 + lam |s|_1 subject to s_t = c_t-g c_{t-1} >= 0

%% inputs:
%   y:  T*1 vector, One dimensional array containing the fluorescence intensities
%withone entry per time-bin.
%   g:  scalar, Parameter of the AR(1) process that models the fluorescence ...
%impulse response.
%   lam:  scalar, sparsity penalty parameter
%   optimize_b: bool, optimize baseline if True
%   optimize_g: integer, number of large, isolated events to consider for
%       optimizing g
%   decimate: int, decimation factor for estimating hyper-parameters faster
%       on decimated data
%   maxIter:  int, maximum number of iterations
%   active_set: npool x 4 matrix, warm stared active sets
%   gmax:  scalar, maximum value of g
%   tau_range: [tau_min, tau_max]
%% outputs
%   c: T*1 vector, the inferred denoised fluorescence signal at each time-bin.
%   s: T*1 vector, discetized deconvolved neural activity (spikes)
%   b: scalar, fluorescence baseline
%   g: scalar, parameter of the AR(1) process
%   active_set: npool x 4 matrix, active sets

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging


%% input arguments
y = reshape(y, [], 1);
T = length(y);

if ~exist('g', 'var') || isempty(g)
    g = estimate_time_constant(y, 1);
end
if ~exist('lam', 'var') || isempty(lam);   lam = 0; end
if ~exist('smin', 'var') || isempty(smin)
    smin = 0;
elseif smin<0
    smin = abs(smin) * GetSn(y);     % use a threshold that is proportional to the noise level
end
if ~exist('optimize_b', 'var') || isempty(optimize_b)
    optimize_b = false;
end
if ~exist('optimize_g', 'var') || isempty(optimize_g)
    optimize_g = 0;
end
if ~exist('decimate', 'var') || isempty(decimate)
    decimate = 1;
else
    decimate = max(1, round(decimate));
end
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = 10;
end
if ~exist('tau_range', 'var') || isempty(tau_range)
   g_range = [0, 1];
else
   g_range = exp(-1./tau_range); 
   g = min(max(g, g_range(1)), g_range(2)); 
end

% change parameters due to downsampling
if decimate>1
    decimate = 1;  %#ok<NASGU>
    disp('to be done');
    %     fluo = y;
    %     y = resample(y, 1, decimate);
    %     g = g^decimate;
    %     thresh = thresh / decimate / decimate;
    %     T = length(y);
end

%% optimize parameters
if ~optimize_b   %% don't optimize the baseline b
    %% initialization
    b = 0;
    [solution, spks, active_set] = oasisAR1(y, g, lam, smin);
    
    %% iteratively update parameters g
    if  optimize_g     % update g
        [solution, active_set, g, spks] = update_g(y, active_set,lam, smin, g_range);
    end
else
    %% initialization
    b = quantile(y, 0.15);
    [solution, spks, active_set] = oasisAR1(y-b, g, lam, smin);
    
    %% iteratively update parameters g and b
    for m=1:maxIter
        b = mean(y-solution);
        %%%%%%%%%% this is a check on the active set %%%%%%%%%%
        if size(active_set, 1) == 0
            break;
        end
        if  optimize_g     % update g
            g0 = g;
            if g>gmax  % spike counts are too small. stop
                g = estimate_time_constant(y, 1);
                [solution, spks, active_set] = oasisAR1(y-b, g, lam, smin);
                break;
            end
            [solution, active_set, g, spks] = update_g(y-b, active_set,lam, smin, g_range);
            if abs(g-g0)/g0 < 1e-3  % g is converged
                optimize_g = false;
            end

        else
            break;
        end
    end
end
c = solution;
s = spks;
end