% we implement thresholded FISTA for regulrizing Ax=b
% via minimization of J(x) = 2*tau*norm(xn,1) + norm(A*xn-b,2)^2;
% we have A,x,b; norm(A)<1
% 
% here's how we do thresholded FISTA:
% 
% let T(x) = S_tau( x + A'*b - A'*A*x )
% next,
% pick some x0 and iterate:
% x_{n+1}=T(x_n+\frac{t_{n}-1}{t_{n+1}} (x^{(n)}-x^{(n-1)}))
% where {t_n} is a sequence of numbers, generated by:
% t_{n+1}=\frac{1+\sqrt{1+4(t_{n})^2}}{2}$ and t_1=1
% see Beck and Teboulle paper for details
function [w_sol,num_iters_taken,tElapsed] = thr_fista_with_wavelets_solve(A,b,Mtb,alpha,tau,w0,maxiters) 

tStart = tic;
m = size(A,1);
n = size(A,2);

% generate sequence of tn's needed for thresholded FISTA
tns = repmat(0,maxiters+10,1);
tns(1) = 1;
for i=2:length(tns)
    tns(i) = (1 + sqrt(1 + 4*tns(i-1)^2))/2;
end

num_iters_taken = 0;
wn = w0;
wn1 = w0;


fprintf('in thr_fista_solve, starting loop..\n');
fprintf('maxiters = %d, tau = %f\n', maxiters, tau);
for i=1:maxiters
            fprintf('doing FISTA with haar iter %d of %d with tau = %f\n', i, maxiters, tau);

            num_iters_taken = i;

            %x_{n+1}=T(x_n+\frac{t_{n}-1}{t_{n+1}} (x^{(n)}-x^{(n-1)}));
            if i==1
                pn = wn;
            else
                pn = wn + ((tns(i)-1)/tns(i+1))*(wn - wn1);
            end
            
            wn1 = wn; % record w_{n-1}

            vn = (1/alpha)*A*haar_1d_inverse(pn);
            zn = (1/alpha)*haar_1d(A'*vn);
			% this is the ISTA step
            wn = softThreshold(pn + Mtb - zn, tau);
end
fprintf('done with thr_fista_solve.. num_iters_taken = %d\n', num_iters_taken); 

% record solution
w_sol = wn;

% record time
tElapsed = toc(tStart);

