function modeList = modes(dens,start)%% modes = modes(kde [,init]) -- Find modes of a KDE via fixed point iter. scheme%   options:%    init -- initial locations for search (default is kde's kernel centers)%    %% NOTE: this process is not guaranteed to find all modes; while it stands an%   excellent chance for Gaussian mixtures, KDEs consisting of Ep. or Lap. kernels%   have discontinous derivatives, leading to quite "jagged" distributions, and%   may have many more modes than kernel centers. See M. Carreira-Perpinan's webpage,%   http://www.cs.toronto.edu/~miguel/research/GMmodes.html, for an excellent%   discussion.%% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txtif (nargin==1) start = getPoints(dens); end;tol=1e-4; 					% set tolerance values etc. max_it=1000; minDistance = 1e-2;pts = getPoints(dens); wts = getWeights(dens); start = start;Ndim = size(pts,1); Npts = size(pts,2); Nloc = size(start,2);BW = getBW(dens,1:Npts); bwMin = min(BW,[],2);modeList = []; vals = [];if (dens.type == 0)  logBW = log(BW);                              % Cache log bandwidths for efficiencyendfor m=1:Nloc					% From each location given:  x = start(:,m);				%   rename for convenience  xTmp = x+inf;  iter = 1;  % FIXED POINT ITERATION TO FIND A MODE:  while (tol < dist(x,xTmp,bwMin) && iter < max_it)	% Iterate until convergence:    diff = pts - repmat(x,[1,Npts]);		%   get distance from kernel centers     xTmp = x;					%   and compute the update:    if (dens.type == 0)				% GAUSSIAN      K  = exp(sum(-.5*(diff./BW).^2-logBW,1)); %   compute kernel eval'n (missing 2pi)      %K = prod(exp(-.5*(diff./BW).^2)./BW,1);	%   (slower kernel eval'n)      px = wts * K';				%   compute kde eval  ( "" )      x  = pts * (wts .* K)' ./ px;		%   compute recursion    elseif (dens.type == 1)			% EPANETCHNIKOV (Mean-shift like update)      K  = max(1-(diff./BW).^2,0)./BW;		%   compute kernel eval'n      px = wts * prod(K,1)';			%   compute kde eval  ( "" )      for d=1:Ndim,				%   compute update in each dimension        Kd = wts .* prod(K([1:d-1,d+1:Ndim],:),1) .* (K~=0);        Kd = Kd ./ sum(Kd);			%        x(d) = pts(d,:) * Kd';			%      end;    else error('Sorry; KDE type not implemented');    end;        iter = iter + 1;				% increment and continue   end  H = llHess(dens,x);				% Compute &   if max(eig(H)) < 0				% Check the Hessian: if it's    modeList = [modeList,x]; vals = [vals,px];	%  neg. def. it's a mode; save it.  end;end[tmp order] = sort(-vals);			% Sort by descending likelihdmodeList = modeList(:,order);				%m = 1;						% Remove any redundant modes:while (m < size(modeList,2))			%  start with "best" modes and  ind = [m+1:size(modeList,2)];			%  work downwards:  d = dist( modeList(:,ind), modeList(:,m+0*ind), bwMin(:,1+0*ind));	  ok  = find(d > minDistance);			%  remove any within minDistance  modeList = modeList(:,[1:m,m+ok]); 		%  of a better mode   m = m+1;end;function d=dist(x,y,bw)  d = sqrt( sum( ((x-y)./bw).^2 ,1));