function [x, Info] = lyngby_opt_gradient(func, x0, varargin) 

% lyngby_opt_gradient   - Optimization with Gradient minimization
%
%       function [x, Info] = lyngby_opt_gradient(func, x0)
% 
%       Input:    func       Function handle (@function)
%                 x0         Vector with parameters to be optimized
%
%       Property: MaxIterations [ {500} ] Positive integer
%                 Parameters    Additional parameters
%
%       Output:   x          Optimized parameters
%                 Info       Structure with additional information
%
%       Optimization with gradient descent. The step size is adapted
%       depended.
%
%       The input argument 'func' should contain a function handle
%       constructed with the at-operator '@'. The function should
%       return the function value and the gradient.
%
%       Examples: 
%         function [f, df] = sin2(x)
%           f = sin(x);
%           df = cos(x);
%
%         function [f, df] = besseljn(x, n)
%           h = abs(x) * 10^6 * eps;
%           f = besselj(x, n);
%           df = mean(gradient(besselj([x-h x x+h], n), h));
%
%         lyngby_opt_gradient(@sin2, 4)
%         lyngby_opt_gradient(@besseljn, 2.8, 'parameters', {8})
%
%       See also LYNGBY, LYNGBY_OPT_CG. 
% 
% $Id: lyngby_opt_gradient.m,v 1.1 2004/03/08 21:42:21 fnielsen Exp $


    % Check arguments
    if nargin < 2
      error('Wrong number of arguments.');
    end
    

    % Default properties
    info = 0;
    maxIterations = 500;
    maxLineSearch = 20;
    minFuncChange = 1e-6;
    parameters = [];


    % Parse properties
    n = 1;
    while n < nargin-1
      arg = lower(varargin{n});

      if strcmp(arg, 'parameters')
	n = n + 1;
	parameters = varargin{n};

      elseif strcmp(arg, 'info')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg)
	  info = arg;
	else
	  error('Wrong type to ''info'' property');
	end

      elseif strcmp(arg, 'maxiterations')
	n = n + 1;
	arg = varargin{n};
	if length(arg) & isreal(arg) & arg > 0 
	  maxIterations = arg;
	else
	  error('Wrong type to ''MaxIterations'' property');
	end
	  
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    argstr = 'feval(func, X';
    for n = 1:length(parameters)
      argstr = [ argstr sprintf(', parameters{%d}', n) ];
    end
    argstr = [ argstr ');' ];

    
    step = 0.1;

    X = x0;
    [fXold GXold] = eval(argstr);
    Xold = X;
    
    for iteration = 1:maxIterations
      X = Xold - step * GXold;
      [fXnew GXnew] = eval(argstr);
      if fXnew < fXold
	step = step * 1.5;
	Xold = X;
	fXold = fXnew;
      else
	step = step / 2;
	X = Xold;
      end
      if info > 5 
	disp(sprintf('Iteration=%03d; f(x)=%f; step=%5e', iteration, ...
	    fXold, step))
      end
      
    end
    
    x = X;



