function output = TVS_FCSAd(input,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% minimize alpha*||Q*y|| + beta*||Phi*y||_1 + 0.5*||Ry-b||_2^2
%%% Here y=Phi'*x  Phi:DWT,  Phi': IDWT  Q: DCT
%%% Sep. 16, 2009, Created by xxxxxxxx in xxxxxxxxx University 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parsin.MAXITER=1; parsin.tv='l1'; l=input.l; u=input.u;
if((l==-Inf)&&(u==Inf))
    project=@(x)x;
elseif (isfinite(l)&&(u==Inf))
    project=@(x)(((l<x).*x)+(l*(x<=l)));
elseif (isfinite(u)&&(l==-Inf))
     project=@(x)(((x<u).*x)+((x>=u)*u));
elseif ((isfinite(u)&&isfinite(l))&&(l<u))
    project=@(x)(((l<x)&(x<u)).*x)+((x>=u)*u)+(l*(x<=l));
else
    error('lower and upper bound l,u should satisfy l<u');
end
t00 = cputime; t0=t00;
n1=input.n1; n2=input.n2; N=n1*n2;
A=input.A; Phi=input.Phi;b=input.b; Atb=A'*b;
TV_iso = false; TV_eps = 0;
y=zeros(input.n1*input.n2,1);
f_prev=inf;f=inf; 
tnew=1; yr=zeros(size(y));
oitr = 1;  Lx=input.L; tnew=1;
alpha = input.alpha; beta = input.beta; 

iterno=1; 
output.funv(iterno)=0.5*sum(b(:).*b(:));
output.xrate(iterno)=1;
output.snr(iterno)=0;
output.xtime(iterno)=0;

KKK=SS1D_InitK(n1*n2);

% iterno=0;

for itr = 1:input.maxitr  % total iter counter        
    iterno=iterno+1; 
%     iitr = iitr+1;
    told=tnew;
    yp=y;
    yg=yr-(A'*(A*yr)-Atb)/Lx;   
%     yg1=reshape(yg,[n1,n2]);
%     if (itr==1)
%         [y, P]=denoise_TV_One(yg, 2*alpha/Lx,-inf,inf,[],parsin);
%     else
%         [y, P]=denoise_TV_One(yg, 2*alpha/Lx,-inf,inf,P,parsin);
%     end
    tmp=dct(yg);
    tmp = shrink1D(tmp,zeros(N,1),2*alpha.*KKK/Lx);
    y=idct(tmp);
%     y=reshape(y, [n1*n2,1]);
    xg=Phi*yg;
    x = shrink1D(xg,zeros(N,1),2*beta/Lx);
    z=0.5*(y+Phi'*x); z=project(z);
    
    [f, f_prev]=get_f(A*z, z, alpha, beta, f);
    
    if f>f_prev  y=yp; f=f_prev; else y=z; end  
    tnew=(1+sqrt(1+4*told^2))/2; 
%     tnew=1;
    yr=y+(told/tnew)*(z-y)+((told-1)/tnew)*(y-yp);    
%     if itr>30 tnew=1; end
%     if (told-1)/tnew>0.5 tnew=1; end
    
%     output.rel(iterno)=norm(y-yp, 'fro')/norm(yp, 'fro');
    output.funv(iterno)=f;
    output.xrate(iterno)=norm(x-input.xs, 'fro')/norm(input.xs, 'fro');
    output.snr(iterno)=snr(y, input.f);
    output.xtime(iterno)=cputime-t00; 
%     if iterno>2 if output.rel(iterno)>output.rel(iterno-1); tnew=1; end  
%     end
    
     if itr>=input.no  
        break
     end
end

output.y=y; output.x=z;
return;

     function [f, f_prev]=get_f(Ay, y, alpha, beta, f)
        v1 = alpha*TV(y);
        v2 = beta*norm(Phi*y,1);
        r3 = Ay - b; 
        v3 = 0.5*(norm(r3(:)).^2);
        f_prev = f;
        f = v1 + v2 + v3;
     end
   
    function p=Grad1(u)
        % backward finite difference along dim 1
        u = reshape(u,n1,n2);
        p = [u(1,:);diff(u,1,1);];
        p = p(:);
    end

    function q=Grad2(u)
        % backward finite difference along dim 2
        u = reshape(u,n1,n2);
        q = [u(:,1) diff(u,1,2)];
        q = q(:);
    end

    function TVx=TV(x)
        % Total Variation norm of x, x is a n by n matrix
        grad_x = [Grad1(x) Grad2(x)];
        if TV_iso
            pt_sqsum = sum(grad_x.*grad_x,2);
            if TV_eps == 0; TVx = sum(sqrt(pt_sqsum)); else TVx = sum(sqrt(pt_sqsum+TV_eps)); end
        else
            TVx = norm(grad_x(:),1);
        end
    end

    function s1 = shrink1D(p,q,s)
        % shrink p toward q by up to s; p and q are vectors, s is scalar
        tmp = p - q; s1 = q + sign(tmp).*max(0,abs(tmp)-s);
    end
end