function [x, u, z, relative_delta, h, v] = tvd_gme(y, lam, K, plot_txt)

% x = tvd_gme(y, lam, K)
%
% GME-TV denoising
% TV denoising using the generalized Moreau envelope
%
% INPUT
%   y : data
%   lam : regularization parameter, lam > 0
%   K : high-pass filter parameter (K non-negative integer)
%               (K = 0 gives scalar MC penalty)
%
% OUTPUT
%   x : denoised data

% To monitor and verify convergence of the algorithm, 
% use the syntax:
% x = tvd_gme(y, lam, K, 'plots')
% As the algorithm converges, the relative error should go to 
% zero, and the points in the scatter plot should 
% go to the graph of the sign function.

% Ivan Selesnick, 2018

% Reference: 
% 'Non-convex total variation regularization for convex denoising of signals'
% I. Selesnick, A. Lanza, S. Morigi, and F. Sgallari
% J. Math. Imaging and Vision, 2020. DOI: 10.1007/s10851-019-00937-5

% This implementation uses ISTA to solve the subproblem
% (optimization of variable v).

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

STOP_TOL = 1e-6;
Max_Iters = 20;

N = length(y);
CS = @(x) cumsum(x(1:end-1));  % cumulative sum

relative_delta = inf;

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

hh = ones(1,K)/K;
h = [zeros(1,K-1) 1 zeros(1,K-1)] - conv(hh,hh);   % length 2K-1

g = deconv(h, [1 -1]);
G = sparse_convmtx(g, N-1-(length(g)-1));
% G = sparse( convmtx(g, N-1-(length(g)-1)) );

c = (1/sqrt(lam)) * g;
C = (1/sqrt(lam)) * G;

% Initialization
x = tvd_L1(y, lam);
Dx = D * x;
x_old = x;
Dx_old = Dx;

% ISTA parameters
v = zeros(N-1, 1);
Cf = freqz(c);
rho = max(abs(Cf))^2;  % rho is max eigenvalue of frequency response
alpha = 1.6 / rho;     % alpha < 2/rho
Nit_ista = 10;
soft = @(t, T) max(t - T, 0) + min(t + T, 0);

iter = 0;

while (relative_delta > STOP_TOL) && (iter < Max_Iters)    
    iter = iter + 1;
    
    % minimize_{v} 0.5 ||C(Dx-v)||_2^2 + ||v||_1
    % i.e., 0.5*sum(abs(C*(Dx-v)).^2) + sum(abs(v))
    % ISTA:
    for k = 1:Nit_ista
        v = soft(v + alpha*(C'*(C*(Dx - v))), alpha);
    end
    
    z = D'*(C'*(C*(Dx - v)));
    x = tvd_L1(y + lam * z, lam);
    Dx = D * x;
    
    relative_delta = norm(D*(x - x_old), 2) / norm(Dx_old, 2);
    x_old = x;
    Dx_old = Dx;
    
    % Optimality scatter plot
    if SHOW_PLOTS
        % Optimality scatter plot
        u = CS(x - y)/lam + C'*(C*(Dx - v));
        
        figure(3)
        clf
        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, 'bo')
        xlim([-M M])
        ylim([-2 2])
        title(sprintf('Optimality scatter plot [Iteration %.2d]', iter))
        drawnow
        
        fprintf('Iteration %.2d, relative error = %.10f\n', iter, relative_delta)
    end
end

u = CS(x - y)/lam + C'*(C*(Dx - v));

end
