function [rec,histout] = NP(data,fun_params,params)
%NP   Newton Projection method. 
%   NP uses the NP method for the image restoration problem
%                    min   J0(X)+ETA*JR(X)
%                    s.t.  X>=0     
%   where J0(X)=0.5*||A*X-Y||^2 is the least squares term (Y is the noisy 
%   and blurred image and A is the blurring operator) and JR(X)=0.5*||X||^2
%   is the Tikhonov regularization functional. 
% 
%   [REC]=NP(DATA,PARAMS) returns the restored image REC. DATA is a 
%   structure array with the problem data. PARAMS is a structure 
%   with the numerical parameters.
%
%   [REC,HISTOUT]=NP(DATA,PARAMS) also returns an array HISTOUT storing 
%   numerical performance of NP:
%       HISTOUT = [ERR F_OBJ G_PROJ TOT_FFT]
%   where ERR is a column vector with the relative error values at the NP
%   iterates; F_OBJ is a column vector with the objective function values; 
%   G_PROJ is a column vector with the projected gradient norm values;
%   TOT_FFT is a column vector with the number of performed FFTs.
%
%   See also NCGP, QNP.

% Get data 
obj = data.obj; % exact image
gn = data.gn; % blurred and noisy image
TF = data.TF; % fft of the psf image
rec0 = data.rec0; % initial guess

% Set up objective function
eta = fun_params.eta; % regularization parameter

% Get stopping criteria parameters
max_fft = params.max_fft; % maximum number of FFTs
Kmax = params.kmax; % maximum number of iterations
tol_gp = params.tol_gp; % projected gradient tolerance
tol_f = params.tol_f; % objective function tolerance
epsilon = 1e-15; % constant for active set

obj_sum = sum(sum(obj.^2));
S = abs(TF).^2;

% Initialize 
rec = rec0;

% Gradient
resid = real(ifft2(TF.*fft2(rec)))-gn;              
g = -real(ifft2(conj(TF).*fft2(resid)))-eta*rec; % -gradient 

% Projected Gradient
ActiveSet = (rec == 0); InactiveSet = (1-ActiveSet);
gp = -g; % gradient
gp = gp.*(InactiveSet+ActiveSet.*(gp < 0)); % projected gradient         

% Numerical values at the initial iterate
err = sqrt(real(sum(sum((real(rec)-obj).^2)))/obj_sum); % relative error
f_obj = 0.5*resid(:)'*resid(:)+eta*0.5*rec(:)'*rec(:); % objective function
g_proj = norm(gp,'fro'); % projected gradient norm
Tot_FFT = 4; % Total FFTs

stop_flag = 1; k = 0; 
tol_gp = tol_gp*g_proj;
histout = zeros(Kmax,4);
histout(1,:) = [err f_obj g_proj Tot_FFT];

%  NP method
while stop_flag
    
    f0_obj = f_obj;
    
    % Evaluate "active set"
    wk = norm(rec-max(0,rec+g),'fro'); epsilonk = min([epsilon;wk]);
    Ik = (rec<=epsilonk & g<0);
    
    % Compute search direction
    gIk = g; gIk(Ik) = 0; 
    d = real(ifft2(fft2(gIk)./(S+eta))); Tot_FFT = Tot_FFT+2;  
    d(Ik) = g(Ik);

    % line search
    [rec,resid,f_obj,void,cft] = ...
        armijo_NP(rec,d,gn,-g,f_obj,eta,TF,epsilonk); 
    if void == 1, disp(' Armijo failure, too many reductions '); break; end
    Tot_FFT = Tot_FFT+cft; 
    
    % relative error
    err = sqrt(real(sum(sum((real(rec)-obj).^2)))/obj_sum); 
    
    % Update gradient           
    g = -real(ifft2(conj(TF).*fft2(resid)))-eta*rec;   
    Tot_FFT = Tot_FFT+2;

    % Projected Gradient
    ActiveSet = (rec == 0); InactiveSet = (1-ActiveSet);
    gp = -g;                                                
    gp = gp.*(InactiveSet+ActiveSet.*(gp < 0));           
    g_proj = norm(gp,'fro');                                
                                                               
    % Get numerical values
    k = k+1;
    histout(k+1,:) = [err f_obj g_proj Tot_FFT];
    
    % check stopping conditions  
    if k >= Kmax
        fprintf(' Maximum number of iterations reached \n'); 
        stop_flag = 0;
    elseif g_proj <= tol_gp
        fprintf(' Met tolerance on the projected gradient \n'); 
        stop_flag = 0;
    elseif Tot_FFT >= max_fft
        fprintf('Maximum number of iterations reached \n'); 
        stop_flag = 0;
    elseif f0_obj-f_obj <= tol_f*f_obj; 
        fprintf(' *** Met relative tolerance on two successive values \n'); 
        fprintf('      of the objective function *** \n');  
        stop_flag = 0;
    end
    
end
histout = histout(1:k+1,:);

