function [x, cost] = tvd_ncvx_ver2(y, lam, pen, Nit)
% [x, cost] = tvd_ncvx_ver2(y, lam, pen, Nit)
% Total variation denoising with non-convex penalty.
%
% This version uses a direct TVD solver in the MM iteration.
% It calls Condat's direct solver.
%
% INPUT
%   y - noisy signal
%   lam - regularization parameter
%   pen - penalty ('log', 'atan', 'mcp', 'L1')
%
% OUTPUT
%   x - denoised signal
%   cost - cost function history

% Ivan Selesnick, selesi@nyu.edu, 2012
% Revised June 2014
% Revised October 2016 (to include 'mcp' for MC penalty)
% Revised January 2017 (simplify D matrix)

% Algorithm: 
% Set lam1 = lam, lam0 = 0 in equations (4) and (36) of the paper:
% 'Convex fused lasso denoising with non-convex regularization 
% and its use for pulse detection'.
% Ankit Parekh and Ivan Selesnick. 
% IEEE Signal Processing in Medicine and Biology Symposium (SPMB), December 2015.


if isrow(y) 
    col = false;
    y = y';                 % Convert to col vector
else
    col = true;
end
cost = zeros(1, Nit);       % Cost function history
N = length(y);

a = 1 / (4 * lam);

[phi, dphi, wfun, ds] = penalty_funs(pen, a);

I = speye(N);
D = I(2:N, :) - I(1:N-1, :);

% Initialization
x = y;                                                
Dx = D*x;

for k = 1:Nit
    x = tvd((y - lam*D'*(ds(Dx)))', N, lam)';
    Dx = D*x;
    cost(k) = 0.5*sum(abs(x-y).^2) + lam*sum(phi(Dx)); % cost function history
end

if col == false
    x = x';     % Convert to row
end

    