function Y=kdopt2d(A,x0,y0,h,maxi,tol)

% KDOPT2D     Numerical optimization (Steepest descent method) for KDE's.
%          
%             Examples: 
%        
%             (1) load kddata;
%                 kdopt2d([glass(:,5) glass(:,8)],0.5,3,0.5,15,1e-2); 
% 
%             KDOPT2D(A,X0,Y0,H,MAXI,TOL) is an algorithm to find a maximum
%             point of a kernel density estimate defined by data A and smoothing 
%             parameter H by the steepest descent method with (X0,Y0) as 
%             the starting guess.  TOL is the required tolerance.  MAXI is the 
%             maximum number of iterations allowed.  
%
%             Only a single value of H is supported, i.e. smoothing is the same
%             in both directions.
%
%             Christian C. Beardah 1994

if (nargin ~= 6)
  error('Incorrect number of function parameters');
end; % of if statement.

% figure('Name','Numerical Optimization for KDE''s');
% clf;

gradf=zeros(2,1);

[n,d]=size(A);

delta=1e-2;

i=0;

vx(1)=x0;
vy(1)=y0;

fn(1)=kdest2d(vx(1),vy(1),A,'knorm',h);

fprintf('i= %2g    x_i= %-8.5f    y_i= %-8.5f    f(x_i,y_i)= %-8.5f\n',i,vx(1),vy(1),fn(1))

i=1; % i is the iteration counter.

while i <= maxi,

  for j=1:n,
    gradf=gradf+knorm((x0-A(j,1))/h,(y0-A(j,2))/h)*[(x0-A(j,1))/h (y0-A(j,2))/h]';
  end;
  gradf=-gradf/(n*h^(d+1));

  r0=[x0 y0]';

  for j=1:21,
    M(j,:)=(r0+(j-11)*gradf/10)';
  end;  

  [N,I]=max(kdest2d(M(:,1),M(:,2),A,'knorm',h));

  r=M(I,:)'; % Find new iterate.

  vx(i+1)=r(1); % Store iterate.
  vy(i+1)=r(2); 

  fn(i+1)=N;

  fprintf('i= %2g    x_i= %-8.5f    y_i= %-8.5f    f(x_i,y_i)= %-8.5f\n',i,vx(i+1),vy(i+1),fn(i+1))

  if (norm(r-r0) < tol), % Converged?
    Y=r; % Answer.

%    xmax=max(vx);
%    xmin=min(vx);
%    xrange=xmax-xmin;

%    ymax=max(vy);
%    ymin=min(vy);
%    yrange=ymax-ymin;

%    xa=linspace(xmin-xrange/8,xmax+xrange/8,40);
%    ya=linspace(ymin-yrange/8,ymax+yrange/8,40);

%    [x,y]=meshgrid(xa,ya);
%    z=kdest2d(x,y,A,'knorm',h);
%    contour(xa,ya,z,30);
%    hold;
%    plot(vx,vy,'w');
%    plot(vx,vy,'go');

%    pause;

%    close;

    return; % Exit the function and return the answer.
  end; % of if statement.

  i=i+1; % Update the iteration counter.

  x0=r(1);
  y0=r(2);
 
end; % of while statement.

Y=r;
disp('Maximum number of iterations reached - terminating execution.');
