function output = TVS_FCSAm(input,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% minimize alpha*TV(y) + beta*||Phi*y||_1 + 0.5*||Ry-b||_2^2
%%% Here y=Phi'*x  Phi:DWT,  Phi': IDWT
%%% 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);
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;

% iterno=0;
para.flag=1;
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;    
    if (itr==1)
        [y,para]=ComputeNonLocalTV(yg, 2*alpha/Lx,para); ;
%         [y, P]=denoise_TV_One(yg, 2*alpha/Lx,-inf,inf,[],parsin);
    else
        if itr==0 ||  itr==round(input.no/2)
%             if mod(itr,2)==0 
%             if mod(itr,20)==0 || itr==1 || itr==round(input.no/2)
            
            [y,para]=ComputeNonLocalTV(yg, 2*alpha/Lx,para); 
%             para.flag=0;
        else
            if isempty
        [y, P]=denoise_TV_One(yg, 2*alpha/Lx,-inf,inf,P,parsin);
        end
    end
    xg=Phi*yg;
    x = shrink1D(xg,zeros(n1,n2),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);
        tmp=Phi*y;
        v2 = beta*norm(tmp(:),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