function [rec,histout] = QNP(data,fun_params,params)
%QNP   Quasi Newton Projection method.
%   QNP uses the QNP method for the image restoration problem
%                    min   J0(X)+ETA*JR(X)
%                    s.t.  X>=0    
%   where J0(X)=sum(sum(Y.*log(Y./(A*X+BG))+A*X+BG-Y)) is the 
%   Kullback-Leibler divergence of X from Y (Y is the noisy and 
%   blurred image, A is the blurring operator, BG is the background value) 
%   and the regularization term JR(X)=0.5*||X||^2 is the Tikhonov 
%   regularization function.
%
%   [REC]=QNP(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]=QNP(DATA,PARAMS) also returns an array HISTOUT storing 
%   numerical performance of QNP:
%       HISTOUT = [ERR F_OBJ G_PROJ TOT_FFT]
%   where ERR is a column vector with the relative error values at the QNP
%   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, NP.

% Get data 
obj = data.obj; % exact image
gn = data.gn; % blurred and noisy image
bg = data.bg; % background
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
den = real(ifft2(TF.*fft2(rec)))+bg;  
temp = gn./den;
g = real(ifft2(conj(TF).*fft2(temp)))-1-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); 
f_obj = sum(sum(gn.*log(gn./den)+den-gn))+eta*0.5*norm(rec,'fro')^2; 
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];

%  QNP 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
    temp2 = temp./(den);
    gIk = g; gIk(Ik) = 0; 
    dd = mean(mean(temp2)); 
    d = real(ifft2(fft2(gIk/dd)./(S+eta/dd ))); Tot_FFT = Tot_FFT+2; 
    d(Ik) = g(Ik);

    % line search
    [rec,den,f_obj,void,cft] = ...
        armijo_QNP(rec,d,gn,den,-g,f_obj,eta,TF,bg,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
    temp = gn./den;
    g = real(ifft2(conj(TF).*fft2(temp)))-1-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,:);

