function [c, s, b, g, active_set] = foopsi_exp2(y, taus, lam, optimize_b,...
    optimize_kernel, decimate, maxIter, shift, win, tol, mask, smin)

%% Infer the most likely discretized spike train for the model c = conv(s, kernel)
% Solves the sparse non-negative deconvolution problem
%  min 1/2|c-y|^2 + lam |s|_1 subject to c = conv(s, kernel)>= 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
%   shift: integer scalar, number of frames by which to shift window from on run of
%       NNLS to the next, default-100
%   win: integer acalar, window size 
%   tol: scalar, tolerance parameters 
%   maxIter: scalar, maximum number of iterations before termination 
%   mask: T * 1 boolean vector, restrict potential spike times 
%   smin: scalar, minimum spike size 

%% 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
%   taus: 1*2 vector, [tau_d, tau_r]

%% 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('shift', 'var') || isempty(shift)
    shift = 100; 
end
if ~exist('win', 'var') || isempty(win)
    win = 200; 
end

if ~exist('tol', 'var') || isempty(tol)
    tol = 1e-9; 
end
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = []; 
end
if ~exist('mask', 'var') || isempty(mask)
    mask = true(T,1); 
end
if ~exist('smin', 'var') || isempty(smin)
    smin = 0; 
end

if ~exist('kernel', 'var') || isempty(kernel)
    g = estimate_time_constant(y, 2);
    taus = ar2exp(g); 
    kernel = exp2kernel(taus, win); 
end

if ~exist('lam', 'var') || isempty(lam);   lam = 0; 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

% 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] = onnls(y, kernel, lam, shift, win, tol, maxIter, mask, smin); 
     
    %% iteratively update parameters g
    if  optimize_g     % update g
        pars.vals = taus; 
       [pars, s] = update_tau(y, spks, pars, nspk); 
    end
else
    disp('to be done'); 
%     %% initialization
%     b = quantile(y, 0.15);
%     [solution, spks, active_set] = oasisAR1(y-b, g, lam);
%     
%     %% iteratively update parameters g and b
%     for m=1:maxIter
%         b = mean(y-solution);
%         if  optimize_g     % update g
%             g0 = g;
%             [solution, active_set, g, spks] = update_g(y-b, active_set,lam);
%             if abs(g-g0)/g0 < 1e-3  % g is converged
%                 optimize_g = false;
%             end
%         else
%             break;
%         end
%     end
end
c = solution;
s = spks;
end


%update the AR coefficient: g
% function [c, active_set, g, s] = update_g(y, active_set, lam)
% %% inputs:
% %   y:  T*1 vector, One dimensional array containing the fluorescence intensities
% %withone entry per time-bin.
% %   active_set: npools*4 matrix, previous active sets.
% %   lam:  scalar, curret value of sparsity penalty parameter lambda.
% 
% %% outputs
% %   c: T*1 vector
% %   s: T*1 vector, spike train
% %   active_set: npool x 4 matrix, active sets
% %   g: scalar
% 
% %% 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
% 
% %% initialization
% 
% len_active_set = size(active_set, 1);  %number of active sets
% y = reshape(y,[],1);    % fluorescence data
% maxl = max(active_set(:, 4));   % maximum ISI
% c = zeros(size(y));     % the optimal denoised trace
% 
% %% find the optimal g and get the warm started active_set
% g = fminbnd(@rss_g, 0, 1);
% yp = y - lam*(1-g);
% for m=1:len_active_set
%     tmp_h = exp(log(g)*(0:maxl)');   % response kernel
%     tmp_hh = cumsum(h.*h);        % hh(k) = h(1:k)'*h(1:k)
%     li = active_set(m, 4);
%     ti = active_set(m, 3);
%     idx = ti:(ti+li-1);
%     active_set(m,1) = (yp(idx))'*tmp_h(1:li);
%     active_set(m,2) = tmp_hh(li);
% end
% [c,s,active_set] = oasisAR1(y, g, lam, [], active_set);
% 
% %% nested functions
%     function rss = rss_g(g)
%         h = exp(log(g)*(0:maxl)');   % response kernel
%         hh = cumsum(h.*h);        % hh(k) = h(1:k)'*h(1:k)
%         yp = y - lam*(1-g);     % include the penalty term
%         for ii=1:len_active_set
%             li = active_set(ii, 4);
%             ti = active_set(ii, 3);
%             idx = ti:(ti+li-1);
%             tmp_v = max(yp(idx)' * h(1:li) / hh(li), 0);
%             c(idx) = tmp_v*h(1:li);
%         end
%         res = y-c;
%         rss = res'*res;     % residual sum of squares
%     end
% end