function [x, cost] = ogs1(y, K, lam, pen, rho, Nit)

% x = ogs1(y, K, lam)
% Overlapping Group Shrinkage/Thresholding (1D)
%
% Input
%  y   : input signal (1D array)
%  lam : regularization parameter (lam > 0)
%  K   : group size
%
% Output
%  x   : output signal
%
% Optional inputs
%  x = ogs1(y, K, lam, pen)
%  x = ogs1(y, K, lam, pen, rho)
%  x = ogs1(y, K, lam, pen, rho, Nit)
%
% where
%  pen : penalty ('abs', 'log', 'atan', 'rat', 'mcp')  [default = 'atan']
%  rho : normalized non-convexity parameter, 0 <= rho <= 1  [default = 1]
%  Nit : number of iterations  [default = 25]
%
% Note: if pen is 'abs' then rho is ignored.
%
% Use [x, cost] = ogs1(...) for cost function history. 

% Po-Yu Chen and Ivan Selesnick
% NYU-Poly, 2013
% Revised: May 2016

if nargin < 6
    Nit = 25;       % Default value
end
if nargin < 5
    rho = 1;        % Default value
end
if nargin < 4
    pen = 'atan';   % Default value
end

% Error checking
if rho == 0
    pen = 'abs';
elseif (rho < 0) || (rho > 1)
    error('Error: need 0 <= rho <= 1.')
end

a_max = 1/(lam*K);
a = rho * a_max;

switch pen
    case 'abs'
        phi = @(x) abs(x);
        wfun = @(x) abs(x);
        a = 0;
    case 'rat'
        phi = @(x) abs(x)./(1+a*abs(x)/2);
        wfun = @(x) abs(x) .* (1 + a*abs(x)/2).^2;
    case 'log'
        phi = @(x) 1/a * log(1 + a*abs(x));
        wfun = @(x) abs(x) .* (1 + a*abs(x));
    case 'atan'
        phi = @(x) 2/(a*sqrt(3)) *  (atan((1+2*a*abs(x))/sqrt(3)) - pi/6);
        wfun = @(x) abs(x) .* (1 + a*abs(x) + a^2.*abs(x).^2);
    case 'mcp'        
        phi = @(t) (abs(t) - a/2 * t.^2).*(abs(t) <= 1/a) + 1/(2*a)*(abs(t) > 1/a);
        wfun = @(x) abs(x) ./ max(1 - a*abs(x), 0); 
        % Note: divide by zero OK -> divide by inf -> multiply by zero
    otherwise
        error('penalty must be ''abs'', ''log'', ''atan'', ''rat'', ''mcp''')
end


y = y(:);
h = ones(K, 1);
x = y;
cost = zeros(1, Nit);

for it = 1 : Nit
    r = sqrt( conv(abs(x).^2, h, 'full') );
    cost(it) = 0.5 * sum(abs(x-y).^2) + lam * sum(phi(r));
    v = 1 + lam*conv( 1./( wfun(r) ), h, 'valid');
    x = y./v;
end
