error in deconvolution

[mschachter]

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

[zhoupc]

thanks, @mschachter . I updated the submodule oasis_matlab.