function [x,fx,n,xhistory,fhistory] = conjugategradient(fun,xinit,dfun, varargin)
% conjugate gradient minimization
% [x,fx,n,xhistory,fhistory] = conjugategradient(fun,xinit,[gradientfunction],[params])
% fun --- a string or inline; objective function
% xstart --- a starting guess
% dfun -- a string or inline; the gradient of the objective function
%         Use [] to signify that a numerical gradient should be used.
% Remaining arguments are parameters to the function

% Returns: x --- argmin(fun)
%          fx -- min(fun)
%          n -- number of directions searched

fun = fcnchk(fun);
if nargin > 2 
    if ~isempty(dfun)
        dfun = fcnchk(dfun);
    end
else
    dfun = [];
end

% a terminating tolerance
tol = .0000001;
xinit = xinit(:);


% keep track of the number of directions searched
n = 0;
x = xinit;
lastx = xinit+100*tol; %ficticious to initialize
xhistory = []; fhistory = [];

if isempty(dfun)
    gnew = -gradfun(fun, x, varargin{:});
else
    gnew = -feval(dfun, x, varargin{:});
end

%while sum(abs(x - lastx)) > (sum(tol*x) + tol) % keep going
 while any(abs(x - lastx) > (tol*x) + tol) % keep going
    if rem(n,length(x)) == 0
        % reinitialize conjugate gradient
        gdir = gnew(:);
        hdir = gdir;
    end

   lastx = x;
   n = n+1;
   
   % my software
   [x, fx,n1d] = findmin1d(fun,x,hdir,tol, varargin{:});
   
   % or, equivalently, using fminbnd
   %dirfuninternal([],fun,x,hdir,varargin{:});  % initialization
   % [alpha,fx] = fminbnd(@dirfuninternal, 0,100);
   % x = x + alpha*hdir;
   
   xhistory = [xhistory; x(:)'];
   fhistory = [fhistory; fx];
   
   if isempty(dfun)
       gnew = -gradfun(fun, x, varargin{:});
   else
       gnew = -feval(dfun, x, varargin{:});
   end
  
   goo = gdir(:)'*gdir(:);
   if goo == 0 
       break;
   end
   % hnew = gnew + hdir*(gnew'*gnew)/(goo); % Fletcher-Reeves
   hnew = gnew + hdir*((gnew+gdir)'*gnew)/(goo); % Polak-Ribiere
   
   hdir = hnew;
   gdir = gnew;
   
   if n > 100*length(x)
      warning('Too many iterations in each direction.');
      return;
   end
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res = dirfuninternal(alpha,varargin)
% dirfun --- turns a multivariate function into a one-variable
% function along a line specified by x and d;
% Initialization: must be initialized for each direction
%  syntax: fundir([], fun, x, d)
% Use
%   syntax: fundir(alpha, [optional params]
% Designed for use with fminbnd

persistent fun;
persistent x; % the position
persistent d; % the direction

if isempty(alpha)
    fun = varargin{1};
    x = varargin{2};
    d = varargin{3};
else
    res = feval(fun, x + alpha.*d, varargin{:});
end
    