How can I solve this problem? "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.375012e-18"

Hello, I'm a Biology PhD student with limited knowledge of MATLAB. I've been trying to use this set of codes I got from a scientific paper that should denoise and deconvolve spikes in neuronal calcium imaging data. I've tried, but failed, within the last 7 months to reach the author. The codes work for some of my samples, but keep giving the warning for others. Below, I have the entire error message and all 4 codes I'm using. I have also attached pictures of what the final result from a complete, error-free run of the codes look like (Picture 1), and from an interrupted run with the warning message (Picture 2). Also, I attached a raw data sample titled 'IK3331_D5_OP50_21.5C_minimum.csv'. I believe the error originates from avec = Mm\b; in Line 261 of the second code below. How can I resolve this error, please?

Error message

>> CDLars_modified_20220913_1
testing sample 1
testing sample 2
testing sample 3
testing sample 4
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =
8.375012e-18. 
> In lars_regression_noise>calcAvec (line 261)
In lars_regression_noise (line 93)
In constrained_foopsi (line 197)
In MCEM_foopsi (line 33)
In CDLars_modified_20220913_1 (line 51)
 
Index exceeds the number of array elements. Index must not exceed 1.

Error in constrained_foopsi (line 198)
         sp = spikes(1:T);

Error in MCEM_foopsi (line 33)
    [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,options);

Error in CDLars_modified_20220913_1 (line 51)
    [c,b,c1,g,sn,sp] = MCEM_foopsi(y,[],[],[],[],options);
 

This first code is titled CDLars_modified_20220913_1

% ピーク高のthrehold
peak_threshold = 0.2; % ピークと見なされるための最小高さ >= 0
join_spacing = 1; % 1つのピークと判断されるための最大間隔 >= 0
% 処理するサンプルの番号
processed_samples = [1:28];%[1:14 16:27]; % sammple 15は特異行列のエラーが出た
% データファイルの名前を設定。
datafile = 'IK3331_D5_OP50_21.5C_minimum.csv';

% 以下の不要な警告を出さないようにする（単なるcsvの形式の問題なので重要ではない）
% '警告: 有効な MATLAB 識別子になるように変数名が変更されました。元の名前は
%  VariableDescriptions プロパティに保存されています。'
warning('off','MATLAB:table:ModifiedAndSavedVarnames');

% csvファイルからデータを読み込む。
dat_tab = readtable(datafile);
sample_names = dat_tab.Properties.VariableNames;
dat = table2array(dat_tab);

allresults = cell(length(processed_samples),1); % 結果記録用
peakcounts = cell(length(processed_samples),1); % 結果記録用

% サンプルごとに順次試す。
for sample_nr = 1:length(processed_samples)
    sample_number = processed_samples(sample_nr);
    disp(['testing sample ' num2str(sample_number)]);
    result = [];
    peakcount = [];
    
    % sample_number番目のサンプルのデータとサンプル名
    y = dat(:,sample_number);
    if any(isnan(y)) && find(isnan(y),1)>1  % NaNは処理できないので除く
      y = y(1:(find(isnan(y), 1 )-1));
    end
    if all(~isnan(y))
    
    sample_name = sample_names{sample_number};
    T = length(y);
    
    % 生データをプロット
    figure('Position', [300 50 700 700])
    subplot(3,1,1)
    plot(y,'b')
    xlabel('timepoint')
    ylabel('R/Ro')
    title([num2str(sample_number) ': ' sample_name], 'Interpreter','none')
    hold on
    
    % CD最適化実行（LARS法） 
    method = 'lars';%'lars'
    options = struct('p',2,'method',method);
    [c,b,c1,g,sn,sp] = MCEM_foopsi(y,[],[],[],[],options);
    
    % スパイクによるカルシウム推定結果をプロット
    plot(c+b, 'm')
    legend('rawdata','fitted')
    hold off
    
    % スパイクピーク高の計算
    tr = g.TAU(end,1); % decay constant
    td = g.TAU(end,2); % decay constant
    peaks = sp* ((tr/td)^(tr/(td-tr)) - (tr/td)^(td/(td-tr))) /sqrt(g.g(1)^2+4*g.g(2));

    % 推定されたスパイクの位置をプロット
    subplot(3,1,2)
    plot([])
    hold on
    %plot(peaks)
    positions = find(peaks>0);
    if ~isempty(positions)
    for i=1:length(positions)
        pos = positions(i);
        plot([pos pos], [0 peaks(pos)], 'm')
    end
    end
    xlim([1 T])
    title('spike estimation')
    xlabel('timepoint')
    ylabel('spike height')
    hold off
        
    % 有効ピーク選択、ピーク融合、ピーク面積計算
    peakarea = peaks;
    peakarea(peakarea < peak_threshold) = 0; % threshold以下は0にしてしまう
    processed_peak = zeros(size(peakarea)); % 処理済みピークをマークするため
    while max(peakarea(~processed_peak))>0 % 未処理のピークがあるか
        t = find(peakarea == max(peakarea(~processed_peak)),1); % うち最大の値の点=t
        if t+1<=T
        for t1 = (t+1):T % 右方向に見ていき、つながっているところはピークの位置1箇所に集約
            if peakarea(t1) > 0 % ピークがつながっていれば
                peakarea(t) = peakarea(t) + peakarea(t1); % 値をピークの位置に移す
                peakarea(t1) = 0;
            else % if 0
              if join_spacing>=1 && any( peakarea(t1:min(T,t1+join_spacing)) > 0 ) % join_spacing以内の幅のギャップであれば
                  t2 = find( peakarea(t1:min(T,t1+join_spacing)) > 0 , 1, 'last' ) + t1-1;
                  peakarea(t) = peakarea(t) + sum(peakarea(t2:t1)); % 値をピークの位置に移してつなげる
                  peakarea(t2:t1) = 0;
                  t1 = t2;
              else
                  break % ピークが途切れた
              end
            end
        end
        end
        if t-1>=1
        for t1 = (t-1):-1:1 % 左方向に見ていき、つながっているところはピークの位置1箇所に集約
            if peakarea(t1)>0 % ピークがつながっていれば
                peakarea(t) = peakarea(t) + peakarea(t1); % 値をピークの位置に移す
                peakarea(t1) = 0;
            else
              if join_spacing>=1 && any( peakarea(max(1,t1-join_spacing):t1) > 0 ) % join_spacing以内の幅のギャップであれば
                  t2 = find( peakarea(max(1,t1-join_spacing):t1) > 0 , 1 ) + max(1,t1-join_spacing)-1;
                  peakarea(t) = peakarea(t) + sum(peakarea(t2:t1)); % 値をピークの位置に移してつなげる
                  peakarea(t2:t1) = 0;
                  t1 = t2;
              else
                  break % ピークが途切れた
              end    
            end
        end
        end
        if t==1 || t1 == 1 % t=0にあるピークはアーティファクトとして除く。
            peakarea(t) = 0;    
        end
        processed_peak(t) = 1;
    end
    peakposition = find(peakarea > 0);
    peaksize = peakarea(peakposition);
    peakcount = [peakposition peaksize]; % ピークの位置と面積を列挙

    % 有効スパイクの位置と面積をプロット
    subplot(3,1,3)
    plot([])
    hold on
    %plot(peakarea)
    positions = find(peakarea>0);
    if ~isempty(positions)
    for i=1:length(positions)
        pos = positions(i);
        plot([pos pos], [0 peakarea(pos)], 'm')
    end
    end
    xlim([1 T])
    title('refined spikes')
    xlabel('timepoint')
    ylabel('spike area')
    hold off
    
    % 図をpdfファイルで保存
    fig = gcf;
    fig.PaperPositionMode = 'auto';
    print(fig,['sample' num2str(sample_number) '.pdf'],'-dpdf','-r300','-fillpage')
    
    result = struct('c',c,'b',b,'c1',c1,'g',g,'sn',sn,'sp',sp,'peaks',peaks,'peakarea',peakarea);
        
    end
    
    %結果を保存
    peakcounts{sample_nr} = peakcount;
    allresults{sample_nr} = result;
end

%ピーク位置の結果をエクセルファイルに保存
%if exist('peakpositions.xls','file')
    %delete 'peakpositions.xls'
if  exist('peakpositions.txt','file')
    delete 'peakpositions.txt'
end
for sample_nr = 1:length(processed_samples)
    sample_number = processed_samples(sample_nr);
    if isempty(peakcounts{sample_nr})
      outdata = [{sample_names{sample_number}},{''};{'position'},{'peaksize'}];    
    else
      outdata = [{sample_names{sample_number}},{''};{'position'},{'peaksize'};num2cell(peakcounts{sample_nr})];
    writecell(outdata,['peakpositions' num2str(sample_number) '.txt']); %add by Jang, xlswrite代わりに使った。
    end

    % Excelの'A1'形式で列を指定しないといけないのでその文字を作る
    %col = sample_nr*2-1;
    %chr1 = floor((col-1)/26);
    %chr2 = mod(col-1,26)+1;
    %if chr1 == 0
        %colcode = char(chr2+64);
    %else
        %colcode = [char(chr1+64) char(chr2+64)];
    %end
    %xlswrite('peakpositions.xls',outdata,1,[colcode '1']);
end
    
% すべての変数をmatファイルに保存
save('allresults.mat', 'allresults');

The second code is titled lars_regression_noise

function [Ws, lambdas, W_lam, lam, flag] = lars_regression_noise(Y, X, positive, noise)

% run LARS for regression problems with LASSO penalty, with optional positivity constraints
% Author: Eftychios Pnevmatikakis. Adapted code from Ari Pakman


% Input Parameters:
%   Y:           Y(:,t) is the observed data at time t
%   X:           the regresion problem is Y=X*W + noise
%   maxcomps:    maximum number of active components to allow
%   positive:    a flag to enforce positivity
%   noise:       the noise of the observation equation. if it is not
%                provided as an argument, the noise is computed from the
%                variance at the end point of the algorithm. The noise is
%                used in the computation of the Cp criterion.


% Output Parameters:
%   Ws: weights from each iteration
%   lambdas: lambda values at each iteration
%   Cps: C_p estimates
%   last_break:     last_break(m) == n means that the last break with m non-zero weights is at Ws(:,:,n)


verbose=false;
%verbose=true;

k=1;

T = size(Y,2); % # of time steps
N = size(X,2); % # of compartments

maxcomps = N;

W = zeros(N,k);
active_set = zeros(N,k);
visited_set = zeros(N,k);

lambdas = [];

Ws=zeros(size(W,1),size(W,2),maxcomps);  % Just preallocation. Ws may end with more or less than maxcomp columns


%% 

r = X'*Y(:);         % N-dim vector
M = -X'*X;            % N x N matrix 

%% begin main loop

%fprintf('\n i = ');
i = 1;
flag = 0;
while 1
    if flag == 1;
        W_lam = 0;
        break;
    end
  %  fprintf('%d,',i);
    
    %% calculate new gradient component if necessary

    if i>1 && new && visited_set(new) ==0         

        visited_set(new) =1;    % remember this direction was computed

    end

    %% Compute full gradient of Q 

    dQ = r + M*W;
        
    %% Compute new W
    if i == 1
        if positive
            dQa = dQ;
        else
            dQa = abs(dQ);
        end
        [lambda, new] = max(dQa(:));
    
        if lambda < 0
            disp('All negative directions!')
            break
        end
        
    else

        % calculate vector to travel along 
         
%        disp('calc velocity')
        
        [avec, gamma_plus, gamma_minus] = calcAvec(new, dQ, W, lambda, active_set, M, positive);        
        
        % calculate time of travel and next new direction 
                
        if new==0                               % if we just dropped a direction we don't allow it to emerge 
            if dropped_sign == 1                % with the same sign
                gamma_plus(dropped) = inf;
            else
                gamma_minus(dropped) = inf;
            end 
        end
                   

        gamma_plus(active_set == 1) = inf;       % don't consider active components 
        gamma_plus(gamma_plus <= 0) = inf;       % or components outside the range [0, lambda]
        gamma_plus(gamma_plus> lambda) =inf;
        [gp_min, gp_min_ind] = min(gamma_plus(:));


        if positive
            gm_min = inf;                         % don't consider new directions that would grow negative
        else
            gamma_minus(active_set == 1) = inf;            
            gamma_minus(gamma_minus> lambda) =inf;
            gamma_minus(gamma_minus <= 0) = inf;                
            [gm_min, gm_min_ind] = min(gamma_minus(:));

        end

        [g_min, which] = min([gp_min, gm_min]);

        if g_min == inf;               % if there are no possible new components, try move to the end
            g_min = lambda;            % This happens when all the components are already active or, if positive==1, when there are no new positive directions 
        end
                 
        

        % LARS check  (is g_min*avec too large?)
        gamma_zero = -W(active_set == 1)  ./ avec;
        gamma_zero_full = zeros(N,k);
        gamma_zero_full(active_set == 1) = gamma_zero;
        gamma_zero_full(gamma_zero_full <= 0) = inf;
        [gz_min, gz_min_ind] = min(gamma_zero_full(:));
        
        if gz_min < g_min            
            if verbose
                fprintf('DROPPING active weight: %d.\n', gz_min_ind)
            end
            active_set(gz_min_ind) = 0;
            dropped = gz_min_ind;
            dropped_sign = sign(W(dropped));
            W(gz_min_ind) = 0;
            avec = avec(gamma_zero ~= gz_min);
            g_min = gz_min;
            new = 0;
            
        elseif g_min < lambda
            if  which == 1
                new = gp_min_ind;
                if verbose
                    fprintf('new positive component: %d.\n', new)
                end
            else
                new = gm_min_ind;
                fprintf('new negative component: %d.\n', new)
            end
        end
                               
        W(active_set == 1) = W(active_set == 1) + g_min * avec;
        
        if positive
            if any (W<0)
                min(W);
                flag = 1;
                %error('negative W component');
            end
        end
        
        lambda = lambda - g_min;
    end

    
 

%%  Update weights and lambdas 
        
    lambdas(i) = lambda;
    Ws(:,:,i)=W;
    res = norm(Y-X*W,'fro')^2;
%% Check finishing conditions    
    

    if lambda ==0 || (new && sum(active_set(:)) == maxcomps) || (res < noise)       
        if verbose
            fprintf('end. \n');
        end
        
        break
    end

%%   
    if new        
        active_set(new) = 1;
    end

    
    i = i + 1;
end
%% end main loop 

%% final calculation of mus
if flag == 0
    if i > 1
        Ws= squeeze(Ws(:,:,1:length(lambdas)));
        w_dir = -(Ws(:,i) - Ws(:,i-1))/(lambdas(i)-lambdas(i-1));
        Aw = X*w_dir;
        y_res = Y - X*(Ws(:,i-1) + w_dir*lambdas(i-1));
        ld = roots([norm(Aw)^2,-2*(Aw'*y_res),y_res'*y_res-noise]);
        lam = ld(intersect(find(ld>lambdas(i)),find(ld<lambdas(i-1))));
        if numel(lam) == 0  || any(lam<0) || any(~isreal(lam));
            lam = lambdas(i);
        end
        W_lam = Ws(:,i-1) + w_dir*(lambdas(i-1)-lam(1));
    else
        cvx_begin quiet
            variable W_lam(size(X,2));
            minimize(sum(W_lam));
            subject to
                W_lam >= 0;
                norm(Y-X*W_lam)<= sqrt(noise);
        cvx_end
        lam = 10;
    end
else
    W_lam = 0;
    Ws = 0;
    lambdas = 0; 
    lam = 0;
end
end

%%%%%%%%%% begin auxiliary functions %%%%%%%%%%


%%
function [avec, gamma_plus, gamma_minus] = calcAvec(new, dQ, W, lambda, active_set, M, positive)

[r,c] = find(active_set);
Mm = -M(r,r);


Mm=(Mm + Mm')/2;

% verify that there is no numerical instability 
eigMm = eig(Mm);
if any(eigMm < 0)
    min(eigMm)
    %error('The matrix Mm has negative eigenvalues')  
    flag = 1;
end


b = sign(W);
if new
    b(new) = sign(dQ(new));
end
b = b(active_set == 1);

avec = Mm\b;

if positive 
    if new 
        in = sum(active_set(1:new));
        if avec(in) <0
            new;
            %error('new component of a is negative')
            flag = 1;
        end
    end
end

    

one_vec = ones(size(W));

dQa = zeros(size(W));
for j=1:length(r)
    dQa = dQa + avec(j)*M(:, r(j));
end

gamma_plus = (lambda - dQ)./(one_vec + dQa);
gamma_minus = (lambda + dQ)./(one_vec - dQa);

end

The third code is titled constrained_foopsi

function [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,options)
% spike inference using a constrained foopsi approach:
%      min      sum(sp)
%    c,sp,b,c1
%      subject to: sp >= 0
%                   b >= 0
%                  G*c = sp
%                   c1 >= 0
%           ||y-b-c - c_in|| <= sn*sqrt(T)

%   Variables:
%   y:      raw fluorescence data (vector of length(T))
%   c:      denoised calcium concentration (Tx1 vector)
%   b:      baseline concentration (scalar)
%  c1:      initial concentration (scalar)
%   g:      discrete time constant(s) (scalar or 2x1 vector)
%  sn:      noise standard deviation (scalar)
%  sp:      spike vector (Tx1 vector)

%   USAGE:
%   [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,OPTIONS)
%   The parameters b,cin,g,sn can be given or else are estimated from the data

%   OPTIONS: (stuct for specifying options)
%         p: order for AR model, used when g is not given (default 2)
%    method: methods for performing spike inference
%   available methods: 'dual' uses dual ascent
%                       'cvx' uses the cvx package available from cvxr.com (default)
%                      'lars' uses the least regression algorithm 
%                     'spgl1' uses the spgl1 package available from
%                     math.ucdavis.edu/~mpf/spgl1/  (usually fastest)
%   bas_nonneg:   flag for setting the baseline lower bound. if 1, then b >= 0 else b >= min(y)
%   noise_range:  frequency range over which the noise power is estimated. Default [Fs/4,Fs/2]
%   noise_method: method to average the PSD in order to obtain a robust noise level estimate
%   lags:         number of extra autocovariance lags to be considered when estimating the time constants
%   resparse:     number of times that the solution is resparsened (default 0). Currently available only with methods 'cvx', 'spgl'
%   fudge_factor: scaling constant to reduce bias in the time constant estimation (default 1 - no scaling)

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015 

defoptions.p = 2;
defoptions.method = 'cvx';
defoptions.bas_nonneg = 1;              % nonnegativity option for baseline estimation
defoptions.noise_range = [0.25,0.5];    % frequency range over which to estimate the noise
defoptions.noise_method = 'logmexp';    % method for which to estimate the noise level
defoptions.lags = 5;                    % number of extra lags when computing the AR coefficients
defoptions.resparse = 0;                % number of times to re-sparse solution
defoptions.fudge_factor = 1;            % fudge factor for time constants

if nargin < 6
    options = defoptions;
    if nargin < 5
        sn = [];
        if nargin < 4
            g = [];
            if nargin < 3
                c1 = [];
                if nargin < 2
                    b = [];
                end
            end
        end
    end
end
      
if ~isfield(options,'p');  options.p = defoptions.p;  end
if ~isfield(options,'method'); options.method = defoptions.method; end
if ~isfield(options,'bas_nonneg'); options.bas_nonneg = defoptions.bas_nonneg; end
if ~isfield(options,'noise_range'); options.noise_range = defoptions.noise_range; end
if ~isfield(options,'noise_method'); options.noise_method = defoptions.noise_method; end
if ~isfield(options,'lags'); options.lags = defoptions.lags; end
if ~isfield(options,'resparse'); options.resparse = defoptions.resparse; end
if ~isfield(options,'fudge_factor'); options.fudge_factor = defoptions.fudge_factor; end

method = options.method;    
if isempty(b);
    bas_est = 1;
else
    bas_est = 0;
end
if isempty(c1)
    c1_est = 1;
else
    c1_est = 0;
end

y = y(:);
T = length(y);
y_full = y;
mis_data = isnan(y);
E = speye(T);
E(mis_data,:) = [];

if any(mis_data)
    y_full(mis_data) = interp1(find(~mis_data),y(~mis_data),find(mis_data));
end
    

if isempty(sn)
    sn = GetSn(y_full,options.noise_range,options.noise_method);
end
if isempty(g)

    g = estimate_time_constants(y_full,options.p,sn,options.lags);
    %{
    % Commentout by Y. iino
    while max(abs(roots([1,-g(:)']))>1) && options.p < 3
        warning('No stable AR(%i) model found. Checking for AR(%i) model \n',options.p,options.p+1);
        options.p = options.p + 1;
        g = estimate_time_constants(y,options.p,sn,options.lags);
    end
    %}
    if options.p == 5
        g = 0;
    end
    %fprintf('Stable AR(%i) model found \n',options.p);
    % re-adjust time constant values
    rg = roots([1;-g(:)]);
    if ~isreal(rg); rg = real(rg) + .001*randn(size(rg)); end
    rg(rg>1) = 0.95 + 0.001*randn(size(rg(rg>1)));
    rg(rg<0) = 0.15 + 0.001*randn(size(rg(rg<0)));
    pg = poly(options.fudge_factor*rg);
    g = -pg(2:end);
end
if options.bas_nonneg  % lower bound for baseline
    b_lb = 0;
else
    b_lb = min(y);
end

if strcmpi(method,'dual'); method = 'dual';
elseif strcmpi(method,'cvx'); method = 'cvx';
elseif strcmpi(method,'lars'); method = 'lars';
elseif strcmpi(method,'spgl1'); method = 'spgl1';
else fprintf('Invalid choice of method. Using CVX \n'); method = 'cvx';
end

if strcmpi(method,'dual') && any(mis_data)
    warning('Dual method does not support missing data. Switching to CVX');
    method = 'cvx';
end

if options.resparse > 0 && (strcmpi(method,'dual') || strcmpi(method,'lars'))
    warning('Resparsening is not supported with chosen method. Switching to CVX');
    method = 'cvx';
end

pathCell = regexp(path, pathsep, 'split');
g = g(:);
G = spdiags(ones(T,1)*[-g(end:-1:1)',1],-length(g):0,T,T);
gd = max(roots([1,-g']));  % decay time constant for initial concentration
gd_vec = gd.^((0:T-1)');

switch method
    case 'dual'
         v = G'*ones(T,1);
        thr = sn*sqrt(T-sum(mis_data));
        if bas_est; b = 0; end
        if c1_est; c1 = 0; end
        myfun = @(Ald) lagrangian_temporal_gradient(Ald,thr^2,y(~mis_data)-b-c1*gd_vec(~mis_data),bas_est,c1_est);
        c = [G\max(G*y,0);zeros(bas_est);zeros(c1_est)];
        options_dual = optimset('GradObj','On','Display','Off','Algorithm','interior-point','TolX',1e-8);
        ld_in = 10;
        [ld,~,flag] = fmincon(myfun,ld_in,[],[],[],[],0,[],[],options_dual);        
        if (flag == -2) || (flag == -3)
            warning('Problem seems unbounded or infeasible. Try a different method.');
        end
        if bas_est; b = c(T+bas_est); end
        if c1_est; c1 = c(end); end
        c = c(1:T);
        sp = G*c;
    case 'cvx'
        onPath = ~isempty(which('cvx_begin'));
        if onPath
            c = zeros(T,1+options.resparse);
            sp = zeros(T,1+options.resparse);
            bas = zeros(1+options.resparse,1);
            cin = zeros(1+options.resparse,1);
            w_ = ones(T,1);
            for rep = 1:options.resparse+1
                [c(:,rep),bas(rep),cin(rep)] = cvx_foopsi(y,b,c1,sn,b_lb,g,w_,~mis_data);
                sp(:,rep) = G*c(:,rep);                
                w_ = 1./(max(sp(:,rep),0) + 1e-8);
            end
            sp(sp<1e-6) = 0;
            c = G\sp;
            b = bas;
            c1 = cin;
        else
            error('CVX does not appear to be on the MATLAB path. It can be downloaded from cvxr.com \n');
        end
    case 'lars'
         Ginv = E*[full(G\speye(T)),ones(T,bas_est),gd_vec*ones(1,c1_est)];
         if bas_est; b = 0; end
         if c1_est; c1 = 0; end    
         [~, ~, spikes, ~, ~] = lars_regression_noise(y(~mis_data)-b_lb*bas_est - b - c1*gd_vec(~mis_data), Ginv, 1, sn^2*(T-sum(mis_data)));
         sp = spikes(1:T);
         b = (spikes(T+bas_est)+b_lb)*bas_est + b*(1-bas_est);
         c1 = spikes(end)*c1_est + c1*(1-c1_est);
         c = G\sp;
    case 'spgl1'
        onPath = ~isempty(which('spgl1'));
        if onPath
            Gx = @(x,mode) G_inv_mat(x,mode,T,g,gd_vec,bas_est,c1_est,E);
            c = zeros(T,1+options.resparse);
            sp = zeros(T,1+options.resparse);
            bas = zeros(1+options.resparse,1);
            cin = zeros(1+options.resparse,1);
            w_ = ones(T,1);
            for rep = 1:options.resparse+1
                if bas_est; b = 0; w_ = [w_;1e-10]; end
                if c1_est; c1 = 0; w_ = [w_;1e-10]; end
                options_spgl = spgSetParms('project',@NormL1NN_project ,'primal_norm', @NormL1NN_primal,'dual_norm',@NormL1NN_dual,'verbosity',0,'weights',w_);
                [spikes,r,~,~] = spg_bpdn( Gx, y(~mis_data)-b_lb*bas_est - (1-bas_est)*b-(1-c1_est)*c1*gd_vec(~mis_data), sn*sqrt(T-sum(mis_data)),options_spgl);
                c(:,rep) = G\spikes(1:T); %Gx([spikes(1:T);0],1);                                  %% calcium signal
                bas(rep) = b*(1-bas_est) + bas_est*spikes(T+bas_est)+b_lb*bas_est;       %% baseline
                cin(rep) = c1*(1-c1_est) + c1_est*spikes(end);
                sp(:,rep) = spikes(1:T);                                           %% spiking signal
                w_ = 1./(spikes(1:T)+1e-8);
            end
            b = bas;
            c1 = cin;
            %sn = norm(r)/sqrt(T);
        else
            error('SPGL1 does not appear to be on the MATLAB path. It can be downloaded from math.ucdavis.edu/~mpf/spgl1 \n');
        end
end

    function sn = GetSn(Y,range_ff,method)
        % estimate noise level with a power spectral density method
        L=length(Y);
%         if ~isempty(which('pmtm'))
%             [psd_Y,ff] = pmtm(Y,5/2,1000,1);
%         end
        if ~isempty(which('pwelch'));
            [psd_Y,ff]=pwelch(Y,round(L/8),[],1000,1);
        else
            xdft = fft(Y)';
            xdft = xdft(:,1:round(L/2)+1);
            psd_Y = (1/L) * abs(xdft).^2;
            ff = 0:1/L:1/2;
            psd_Y(2:end-1) = 2*psd_Y(2:end-1);
        end
        ind=ff>range_ff(1);
        ind(ff>range_ff(2))=0;
        switch method
            case 'mean'
                sn=sqrt(mean(psd_Y(ind)/2));
            case 'median'
                sn=sqrt(median(psd_Y(ind)/2));
            case 'logmexp'
                sn = sqrt(exp(mean(log(psd_Y(ind)/2))));
        end
    end

    
    function g = estimate_time_constants(y,p,sn,lags)
        % estimate time constants from the sample autocovariance function
        
        lags = lags + p;
        if ~isempty(which('xcov')) %signal processing toolbox
            xc = xcov(y,lags,'biased');
        else
            ynormed = (y - mean(y));
            xc = nan(lags + 1, 1);
            for k = 0:lags
                xc(k + 1) = ynormed(1 + k:end)' * ynormed(1:end - k);
            end
            xc = [flipud(xc(2:end)); xc] / numel(y);
        end
        xc = xc(:);
        A = toeplitz(xc(lags+(1:lags)),xc(lags+(1:p))) - sn^2*eye(lags,p);
        g = pinv(A)*xc(lags+2:end);            
    end

    function [f,grad] = lagrangian_temporal_gradient(Al,thr,y_raw,bas_flag,c1_flag)
        options_qp = optimset('Display','Off','Algorithm','interior-point-convex');
        H = [speye(T),ones(T,bas_flag),gd_vec*ones(1,c1_flag);...
            ones(bas_flag,T),T*ones(bas_flag),(1-gd^T)/(1-gd)*ones(c1_flag,bas_flag);...
            (gd_vec*ones(1,c1_flag))',(1-gd^T)/(1-gd)*ones(c1_flag,bas_flag),(1-gd^(2*T))/(1-gd^2)*ones(c1_flag,c1_flag)];
        Ay = [y_raw;sum(y_raw)*ones(bas_flag);gd_vec'*y_raw*ones(c1_flag)];
        c = quadprog(2*Al(1)*H,[v;zeros(bas_flag+c1_flag,1)]-2*Al(1)*Ay,[-G,sparse(T,bas_flag+c1_flag);sparse(bas_flag+c1_flag,T),-speye(bas_flag+c1_flag)]...
            ,[sparse(T,1);-b_lb*ones(bas_flag);zeros(c1_flag)],[],[],[],[],c,options_qp);
        f = v'*c(1:T);    
        grad = [sum((c(1:T)-y_raw + c(T+bas_flag)*bas_flag + c(end)*gd_vec*c1_flag).^2)-thr];
        f = f + Al(:)'*grad;
    end

    function b = G_inv_mat(x,mode,NT,gs,gd_vec,bas_flag,c1_flag,Emat)
        if mode == 1
            b = filter(1,[1;-gs(:)],x(1:NT)) + bas_flag*x(NT+bas_flag) + c1_flag*gd_vec*x(end);
            b = Emat*b;
           %b = G\x(1:NT) + x(NT+bas_flag)*bas_flag + x(end)*c1_flag;
        elseif mode == 2
            x = Emat'*x;
            b = [flipud(filter(1,[1;-gs(:)],flipud(x)));ones(bas_flag,1)*sum(x);ones(c1_flag,1)*(gd_vec'*x)];
           %b = [G'\x;ones(bas_flag,1)*sum(x);ones(c1_flag,1)*(gd_vec'*x)] ;
        end
    end

end

The fourth code is titled MCEM_foopsi

function [c,b,c1,g,sn,sp] = MCEM_foopsi(y,b,c1,g,sn,options)
    defoptions.dt = 1;
    defoptions.MaxIter = 10;
    defoptions.MaxInerIter = 50;
    defoptions.TauStd = [0.2,2];
    defoptions.default_g = [0.6,0.9];
    
    if nargin < 6; options.defoptions; end 
    if ~isfield(options,'dt'); options.dt = defoptions.dt; end
    if ~isfield(options,'MaxIter'); options.MaxIter = defoptions.MaxIter; end
    if ~isfield(options,'MaxInerIter'); options.MaxInerIter = defoptions.MaxInerIter; end
    if ~isfield(options,'TauStd'); options.TauStd = defoptions.TauStd; end
    if ~isfield(options,'default_g'); options.default_g = defoptions.default_g; end
    
    dt = options.dt;
    
    if nargin < 5
        sn = [];
        if nargin < 4
            g = [];
            if nargin < 3
                c1 = [];
                if nargin < 2
                    b = [];
                end
            end
        end
    end

    
% initialization    
    
    [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,options);
    T = length(y);
    %p = length(g);
    p = 2; % changed by YIino
    gr_ = sort(roots([1,-g(:)']));
    if p == 1; gr_ = [0,gr_]; end
    if any(gr_<0) || any(~isreal(gr_))
        gr_ = options.default_g;
    end
    tau = -dt./log(gr_);
    tau = tau(1:2); % added by Y. Iino 20180312 because of error in different vector size
    if p == 1        
        gr_(1) = 0;
        G1sp = zeros(T,1);
    else
        %G1sp = make_G_matrix(T,gr_(1))\sp(:);
        % modified by Y. Iino 20180312 because make_G_matrix is missing
        G1 = spdiags(ones(T,1)*[-min(gr_),1],-1:0,T,T);
        G1sp = G1\sp(:);
    end
    tau1_std = max(tau(1)/5,options.TauStd(1));
    tau2_std = min(tau(2)/10,options.TauStd(2));
    tauMoves = [0 0];
    tau_min = 0;
    tau_max = 2*tau(2);    
    gd_vec = max(gr_).^(0:T-1)';
    tau_sam = zeros(options.MaxIter,2);
    for iter = 1:options.MaxIter
        TAU = zeros(options.MaxInerIter,2);
        for iner_iter = 1:options.MaxInerIter
            if p >= 2       % update rise time constant
                logC = -norm(y(:) - c(:) - b - c1*gd_vec(:))^2; 
                tau_ = tau;
                tau_temp = tau_(1)+(tau1_std*randn); 
                while tau_temp >tau(2) || tau_temp<tau_min
                    tau_temp = tau_(1)+(tau1_std*randn);
                end 
                tau_(1) = tau_temp;
                gr_ = exp(dt*(-1./tau_));
                h_ = gr_(2)-gr_(1);
                G1_ = spdiags(ones(T,1)*[-min(gr_),1],-1:0,T,T);
                G1sp_ = G1_\sp(:);
                G2_ = spdiags(ones(T,1)*[-max(gr_),1],-1:0,T,T);
                c_ = (-G1sp_*gr_(1)+G2_\sp(:)*gr_(2))/h_;

                logC_ = -norm(y(:)-c_(:)-b-c1*gd_vec)^2;

                %accept or reject
                prior_ratio = 1;
        %         prior_ratio = gampdf(tau_(2),12,1)/gampdf(tau(2),12,1);
                ratio = exp((logC_-logC)/(2*sn^2))*prior_ratio;
                if rand < ratio %accept
                    tau = tau_;
                    h = h_; G1sp = G1sp_; G2 = G2_; c = c_; %#ok<NASGU> 
                    tauMoves = tauMoves + [1 1];
                else
                    tauMoves = tauMoves + [0 1];
                end                
            end
            %%%%%%%%%%%%%%%%%%%%%%%
            % next update decay time constant
            %%%%%%%%%%%%%%%%%%%%%%%

            %initial logC
            logC = -norm(y(:)-c(:)-b-c1*gd_vec)^2;
            tau_ = tau;
            tau_temp = tau_(2)+(tau2_std*randn);
            while tau_temp>tau_max || tau_temp<tau_(1)
                tau_temp = tau_(2)+(tau2_std*randn);
            end  
            tau_(2) = tau_temp;

            gr_ = exp(dt*(-1./tau_));
            gd_vec_ = max(gr_).^(0:T-1)';
            h_ = gr_(2)-gr_(1);
            %G1_ = spdiags(ones(T,1)*[-min(gr_),1],[-1:0],T,T);
            G2_ = spdiags(ones(T,1)*[-max(gr_),1],-1:0,T,T);  
            c_ = (-G1sp(:)*gr_(1)+G2_\sp(:)*gr_(2))/h_;

            logC_ = -norm(y(:)-c_(:)-b-c1*gd_vec_)^2;

            %accept or reject
            prior_ratio = 1;
    %       prior_ratio = gampdf(tau_(2),12,1)/gampdf(tau(2),12,1);
            ratio = exp((1./(2*sn^2)).*(logC_-logC))*prior_ratio;
            if rand<ratio %accept
                tau = tau_;
                %h = h_; G1 = G1_; G2 = G2_; 
                c = c_; gd_vec = gd_vec_;
                tauMoves = tauMoves + [1 1];
            else
                tauMoves = tauMoves + [0 1];
            end   
            TAU(iner_iter,:) = tau;
        end
        tau_ = mean(TAU);
        tau_sam(iter,:) = tau_;
        gr_ = exp(dt*(-1./tau_));
        %res = y(:)-c_(:)-b-c1*gd_vec;
        %sn   = 1./sqrt(gamrnd(1+T/2,1/(0.1 + sum((res.^2)/2))));
        g = [sum(gr_);-prod(gr_)];
        [c,b,c1,~,sn,sp] = constrained_foopsi(y,[],[],g(1:p),sn,options);
        %disp(tauMoves);
        
    end
    %g.g = g;
    %g.TAU = tau_sam(:,3-p:2);
    % modified by Y. Iino, 20180312 because following error occurred
    % 非構造体配列オブジェクトにフィールドを代入しました。
    g = struct('g',g,'TAU',tau_sam(:,3-p:2));
end

[Torsten]

Who else but the authors of the code can answer your questions ? So try to contact them.

[BintaMaria]

It's a different question if you can't offer any help. You can't think I would skip contacting the author first only to post my question here...

I forgot to include in my description that I've tried several times within the last 8 months to reach the author, all to no avail. I still haven't gotten any responses to my emails or LinkedIn messages.

[Torsten]

But I can assure you that nobody else - except maybe a friend of yours - will dive into this code and search for the possible reason for its failure.

[BintaMaria]

@Torsten speak for yourself and kindly stop commenting on my post if you have nothing helpful to contribute.

[Torsten]

Well, I thought it was helpful to show you the limits of helpfulness.

[Stephen23]

"I believe the error originates from avec = Mm\b; in Line 261 of the second code below."

No, that is just a warning.

While the paying attention to warnings is important, it does not (immediately) stop the code from running.

"How can I resolve this error, please?"

The error you get is caused on this line:

sp = spikes(1:T);

with the stated cause of: "Index exceeds the number of array elements. Index must not exceed 1." So apparently the code's author expected the function lars_regression_noise to return as its third output an array with atleast T elements in it, but in fact T only has one element. Thus the error.

The fact that the function inputs&outputs listed in its help do not match its function signature is not a good sign: apparently someone modified this from someone else's code, but did not bother to update the help to match. Are all of the help and comments wrong, or just part of them? We simply don't know.

In any case, you need to debug lars_regression_noise and find out under what conditions W_lam is scalar or non-scalar. It might be due to the singular matrix, or it might be due to e.g. the data scaling, number of peaks, etc. To start debugging, take a close look at that WHILE loop, its flags, etc. Some likely places for you to investigate:

if lambda < 0
    disp('All negative directions!')
    break
end

and

%% final calculation of mus
if flag == 0
    ..
else
    W_lam = 0;
    ..
end

Those might cause the function to return with a scalar W_lam, thus triggering that error.

[BintaMaria]

Hi Stephen, thank you so much for taking the time assess the problem and giving suggestions. I will do as you have suggested and revert back.

Answer 1

[Matt J]

I believe the error originates from avec = Mm\b

If so, it is because the Mm matrix is singular. The solution would be to make it non-singular. But why it should be non-singular in the first place is something only you or the author of the code would know.

[BintaMaria]

Hi Matt, thanks for your response. If I were to make it non-singular, what would you kindly suggest?

[Matt J]

If you were to make it non-singular, the line avec=Mm\b shouldn't give you a problem anymore.