function [rec,histout] = NCGP(data,fun_params,params)
%NCGP   Newton-Conjugate Gradient Projection method.
%   NCGP uses the NCGP method for the image restoration problem
%                    min   J0(X)+ETA*JR(X)
%                    s.t.  X>=0    
%   where J0 is the data fitting function and JR is the regularization 
%   function. The inner linear system is solved by the Conjugate Gradient 
%   method.
% 
%   [REC]=NCGP(DATA,FUN_PARAMS,PARAMS) returns the restored image REC. 
%   DATA is a structure array with the problem data. FUN_PARAMS is a 
%   structure array with the objective function data. PARAMS is a structure
%   with the numerical parameters.
%
%   [REC,HISTOUT]=NCGP(DATA,FUN_PARAMS,PARAMS) also returns an array
%   HISTOUT storing numerical performance of NCGP:
%       HISTOUT = [ERR F_OBJ G_PROJ TOT_FFT CG_IT]
%   where ERR is a column vector with the relative error values at the NCGP
%   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; 
%   CG_IT is a column vector with the number of CG iterations.
%
%   See also QNP, NP.

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

% Set up objective function
fun_obj = fun_params.obj; % objective function
eta = fun_params.eta; % regularization parameter 
Amult = fun_params.Amult; % function performing the matrix-vector 

% 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
pcg_max = params.pcg_max; % Maximum number of CG iterations
pcg_tol = params.pcg_tol; % CG relative tolerance
precond = params.precond; % if precond=1, CG with diagonal preconditioner
                          % else (precond=0), CG without preconditioner

fprintf('\n-------------------------------------------------\n');
switch fun_obj
    case 'KLTKb_obj'   
        fprintf(' NCGP method for deblurring Poisson data\n');
        fprintf(' with Tikhonov regularization \n');
    case 'KLTVb_obj'   
        fprintf(' NCGP method for deblurring Poisson data\n');
        fprintf(' with Total Variation regularization \n');
    case 'LSTKb_obj'   
        fprintf(' NCGP method for deblurring Gaussian data\n');
        fprintf(' with Tikhonov regularization \n');
    case 'LSTVb_obj'   
        fprintf(' NCGP method for deblurring Gaussian data\n');
        fprintf(' with Total Variation regularization \n');
    case 'KLTVn_obj'   
        fprintf(' NCGP method for denoising Poisson data\n');
        fprintf(' with Total Variation regularization \n');
    case 'LSTVn_obj'   
        fprintf(' NCGP method for denoising Gaussian data\n');
        fprintf(' with Total Variation regularization \n');
end
fprintf('-------------------------------------------------\n\n');

% ----------------------------------------------------------------------
% Initialization
% ----------------------------------------------------------------------
rec = rec0;                                 

% Set up discretization of first derivative operators
[nx,ny] = size(rec); 
if strcmp(fun_obj,'KLTVb_obj') || strcmp(fun_obj,'LSTVb_obj') || ...
        strcmp(fun_obj,'KLTVn_obj') || strcmp(fun_obj,'LSTVn_obj')        
    fun_params.betatv = 0.1;
    fun_params = SetTV(fun_params,nx,ny);
elseif strcmp(fun_obj,'KLTKb_obj') || strcmp(fun_obj,'LSTKb_obj')  
    fun_params.nx = nx; fun_params.ny = ny; 
else
    fprintf('Error in objective function name'); 
end

% Set up parameter for active set
epsilon = 1e-15;  
Tot_FFT = 0;   
obj_sum = sum(sum(obj.^2));

% Compute objective function, its gradient and hessian 
[f_obj,g,Hfit,Hreg,objft] = feval(fun_obj,rec,fun_params,data);
Tot_FFT = Tot_FFT +objft;
g = -g; % -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
g_proj = norm(gp,'fro'); % projected gradient norm

stop_flag = 1; k = 0; cg_it = 0;
tol_gp = tol_gp*g_proj;
histout = zeros(Kmax,5);

% NCGP 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 );
    
    % Right hand side
    gIk = g; gIk(Ik) = 0; 
    
    % Conjugate gradient iterations 
    [d,iters,cgft] = ...
        PCGsolver(gIk,TF,Hfit,Hreg,eta,Amult,pcg_tol,pcg_max,precond); 
    d = reshape(d,nx,ny); 
    cg_it = cg_it+iters; 
    Tot_FFT = Tot_FFT+cgft; 
    d(Ik) = g(Ik);

    % Line search
    [rec,void,cft] = ...
        armijo_NCGP(rec,d,f_obj,-g,data,fun_params,epsilonk); 
    if void == 1, disp(' Armijo failure, too many reductions '); break; end
    Tot_FFT = Tot_FFT +cft;
    
    % Compute objective function, its gradient and hessian 
    [f_obj,g,Hfit,Hreg,objft] = feval(fun_obj,rec,fun_params,data); 
    Tot_FFT = Tot_FFT+objft;
    g = -g; % -gradient
 
    % Projected Gradient
    ActiveSet = (rec == 0); InactiveSet = (1 - ActiveSet);
    gp = -g; % gradient
    gp = gp.*(InactiveSet+ActiveSet.*(gp < 0)); % projected gradient  
    g_proj = norm(gp,'fro'); % projected gradient norm  
         
    % Relative error
    err = sqrt(real(sum(sum((real(rec)-obj).^2)))/obj_sum); 
    
    % Get numerical values
    k = k+1;
    histout(k,:) = [err f_obj g_proj Tot_FFT cg_it];

    % 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 FFTs 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,:);
