function [x, u, cost] = mtvd(y, lam, gamma, Nit, plot_txt)

% x = mtvd(y, lam, gamma, Nit)
%
% Moreau-enhanced Total Variation (TV) Denoising
%
% INPUT
%   y : data
%   lam : regularization parameter, lam > 0
%   gamma : non-convexity parameter, 0 <= gamma < 1
%
% OUTPUT
%   x : denoised data
%
% Suggested: set 0.5 <= gamma <= 0.8

% Reference:
% Total Variation Denoising via the Moreau Envelope
% Ivan Selesnick
% IEEE Signal Processing Letters
% 2017

SHOW_PLOTS = false;
if nargin > 4
    if strcmp(plot_txt, 'plots')
        SHOW_PLOTS = true;
    end
end

N = length(y);

cost = nan(1, Nit);

% Initialization
x = tvd(y, N, lam)';

alpha = gamma / lam;

if gamma > 0
    
    for iter = 1:Nit
        
        z = y + lam * alpha * (x - tvd(x, N, 1/alpha)');
        x = tvd(z, N, lam)';
        
        v = tvd(x, N, 1/alpha)';
        cost(iter) = 0.5 * sum((y-x).^2) + lam * sum(abs(diff(x))) - ...
            lam*sum(abs(diff(v))) - lam*alpha/2 * sum((x-v).^2);
        
        % Plots
        if SHOW_PLOTS
            figure(2)
            clf
            
            % Optimality scatter plot
            u = cumsum( (x - y)/lam + alpha*(tvd(x,N,1/alpha)' - x) );
            u = u(1:end-1);
            
            subplot(2, 1, 1)
            M = 1.1*max(abs(diff(x)));
            if M == 0, M = 1; end
            plot([-M 0 0 M], [-1 -1 1 1], 'r-', diff(x), u, 'ko')
            xlim([-M M])
            ylim([-2 2])
            title(sprintf('Optimality scatter plot [Iteration %.2d]', iter))
            
            subplot(2, 1, 2)
            plot(x)
            xlim([0 N])
            title(sprintf('Denoised signal [Iteration %.2d]', iter))
            drawnow
            
            % Set to true to print plots to a file.
            if false
                orient tall
                if iter == 1
                    print -dpsc mtvd_iters
                else
                    print -dpsc -append mtvd_iters
                end
            end
            
        end
        
    end
    
    u = cumsum( (x - y)/lam + alpha*(tvd(x,N,1/alpha)'-x) );
    u = u(1:end-1);
    
else
    
    u = cumsum( (x - y)/lam );
    u = u(1:end-1);
    
end

