function [x,f]=minimize_conjugate_gradient(x0,fg,optparam)
%
% Jean-Philippe Vert
% June 4, 2006
%
% function [xn,fn]=minimize_conjugate_gradient(x0,fg,optparam)
%
% Minimization of a nonlinear function by conjugate gradient.
% Warning: the line search is backtracking, which is not a good choice in
% thid case.
%
% Input:    x0          initial point
%           fg         name of a matlab function [f,g,H]=fg(x) which returns the
%                       objective function and its gradient and Hessian
%           optparam    optional parameters
%           optparam.tol            stopping tolerance: the algorithm stops when
%                                   ||g(x)||<=tol*min(1,||g(x0)||)
%           optparam.armijo_alpha   alpha parameter for Armijo rule
%           optparam.armijo_beta    beta parameter for Armijo rule
%           optparam.display        display information is =1%
%           optparam.plot           plot the successive iterates (x[1] and
%                                   x[2])
%
% Output:   x      final point after optimization
%           f      objective at the final point

if nargin<3
    optparam=0;
end

if isfield(optparam,'armijo_alpha')
    armijo_alpha=optparam.armijo_alpha;
else
    armijo_alpha=0.1;
end
if isfield(optparam,'armijo_beta')
    armijo_beta=optparam.armijo_beta;
else
    armijo_beta=0.5;
end
if isfield(optparam,'tol')
    tol=optparam.tol;
else
    tol=1e-6;
end
if isfield(optparam,'display')
    display=optparam.display;
else
    display=0;
end
if isfield(optparam,'plot')
    plotit=optparam.plot;
else
    plotit=0;
end
if plotit
    hold on
end

% Initializations
x=x0;
[f,g]=feval(fg,x);
nmg0=norm(g);
nmg=nmg0;
gg=nmg^2;
gam=0;
h=zeros(size(x0));
nh=0;
it=0;
a1=0.001;


% main loop
while (nmg>tol*min(1,nmg0))
    it=it+1;
    if display
        fprintf('it=%3.d  f=%e   ||g||=%e\n',it,f,nmg);
    end
    
    % Previous values
    hpr=h;
    nhpr=nh;
    ggpr=gg;
    gpr=g;
    
    % New direction and scaling
    h = gam*h - g;
    nh=norm(h);
    
    % Check descent
    gh=dot(g,h);
    if gh >= - a1 * nmg * nh
        % Not descent enough: take steepest descent
        fprintf('(taking steepest descent)')
        h = -g;
        nh = nmg;
    end

    % Compute the new point with Armijo rule
    [t,xnew,fnew] = armijo(x,h,gh,fg,1,f,armijo_alpha,armijo_beta);
    if plotit
        plot([x(1),xnew(1)],[x(2),xnew(2)],'-')
    end
    x = xnew;
    
    % Compute the new gradient and update x and gamma
    [f,g] = feval(fg,x);
    nmg = norm(g);
    gg = nmg^2;
    gam = dot(g-gpr,g)/ggpr; % Polak-Ribiere rule
    
end
it=it+1;
if display
    fprintf('it=%3.d  f=%e   ||g||=%e\n',it,f,nmg);
    fprintf('Successful termination!')
end
