% created by Nils Hase 2016 (nilshase@math.uni-bremen.de)

% Fast Iterative Shrinkage Thresholding coupled with the Morozov parameter
% choice rule for problems of the type
%
% x^* = arg min_{x} 1/2*||Ax-y||_2^2 + alpha/2*||x||_2^2
% s.t. x >= 0

% Almost identical code to FISTA_L2, but with projection step,
% differs from FISTA_L1_pos by the shrinkage step

function [x_opt,alpha_opt] = FISTA_L2_pos(A,y,x_0,alpha,plotter)

delta = 1;

[m,n] = size(A);
lambda = 1/norm(A,'fro')^2;
max_iter = 200;     % 200-300 iterations are typically enough
tol = 10^-6;
tau = 1;

discrepancy = zeros(size(alpha));

x_opt = x_0;
for k = 1:length(alpha)
    % Resetting it to x_0 is advantageous (this has been tested)
    x = x_0;
    x_old = x_0;
    x_old_old = x_0;
    
    for kk = 1:max_iter
        x_old_old = x_old;
        x_old = x;
        
        % FISTA-update
        xx = x_old + (kk-2)/(kk+1)*(x_old - x_old_old);
        x = xx - lambda*(A'*(A*xx - y));
        
        % Shrinkage-step (L2)
        x = x/(1+lambda*alpha(k));

        % Projection step
        x(x<0) = 0;
        
        if norm(x_old - x) < tol
            break;
        end
    end

    discrepancy(k) = norm(A*x - y)/sqrt(length(y));
      
    x_opt = x;
    alpha_opt = alpha(k);
    % Morozov's discrepancy principle
    if discrepancy(k) < tau*delta
        fprintf('   Discrepancy principle found optimal alpha!\n')
        break;
    end
end

if plotter == 1
    figure();
    plot(log10(alpha),discrepancy)
    hold on;
    plot([log10(alpha(1)),log10(alpha(end))],tau*delta*[1,1],'r--')
    title('Discrepancy curve FISTA L2 pos')
    hold off;
end

end
