%PARZENML Optimum smoothing parameter in Parzen density estimation.
% 
%   H = PARZENML(A)
% 
% INPUT	
%   A    Input dataset
%
% OUTPUT
%   H    Scalar smoothing parameter (in case of crisp labels)
%        Vector with smoothing parameters (in case of soft labels)
%
% DESCRIPTION
% Maximum likelihood estimation for the smoothing parameter H in the 
% Parzen denstity estimation of the data in A. A leave-one out 
% maximum likelihood estimation is used. 
%
% The dataset A can either be crisp or soft labeled. In case of crisp
% labeling the class information is not used and a single smoothing 
% parameter is estimated. In case of soft labels a smoothing parameter
% for every class is estimated and objects are weighted in relation to
% their class weigthts (soft label value). 
% It may be profitable to scale the data before calling it. eg. 
% WS = SCALEM(A,'variance'); A = A*WS.
% 
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, MAPPINGS, SCALEM, SELDAT, PARZENM, PARZENDC, PRPROGRESS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: parzenml.m,v 1.11 2010/03/25 15:39:46 duin Exp $

function h = parzenml(A,fid)

		
	if nargin < 2, fid = []; end

	if isdouble(A), A = prdataset(A); end
	
	A = testdatasize(A);
	A = testdatasize(A,'objects');

	if islabtype(A,'crisp')
		h = parzenmlc(A,fid);
	elseif islabtype(A,'soft')
		h = parzenmls(A,fid);
	else
		error('Label type should be either ''crisp'' or ''soft''')
	end
	
	return
	
function h = parzenmlc(A,fid) %crisp version

	[m,k] = size(A);
	DD= distm(+A) + diag(1e70*ones(1,m));
	E = min(DD);
	
	h1 = sqrt(max(E));    % initial estimate of h
	F1 = derlc(DD,E,h1,k); % derivative

	prprogress(fid,'parzenml:\n');
	prprogress(fid,' %6.4f   %6.3e\n',h1,F1);
	if abs(F1) < 1e-70 
		h = h1;
		prwarning(4,'jump out\n');
		return;
	end
	
	a1 = (F1+m*k)*h1*h1;
	h2 = sqrt(a1/(m*k));  % second guess
	F2 = derlc(DD,E,h2,k); % derivative

	prprogress(fid,' %6.4f   %6.3e\n',h2,F2);
	if (abs(F2) < 1e-70) | (abs(1e0-h1/h2) < 1e-6) 
		h = h2;
		prwarning(4,'jump out\n');
		return
	end
	
	% find zero-point of derivative to optimize h^2
	% stop if improvement is small, or h does not change significantly
	
	alf = 1;
	prwaitbar(100,'parzenml: Optimizing smoothing parameter',m > 100)
	iter = 0;
	while abs(1e0-F2/F1) > 1e-4 & abs(1e0-h2/h1) > 1e-3 & abs(F2) > 1e-70
		iter = iter+1;
		h3 = (h1*h1*h2*h2)*(F2-F1)/(F2*h2*h2-F1*h1*h1);
		if h3 < 0 % this should not happen
			h3 = sqrt((F2+m*k)*h2*h2/(m*k));
		else
			h3 = sqrt(h3);
		end
		prwaitbar(100,100-100*exp(-iter/10));
		h3 = h2 +alf*(h3-h2);
		F3 = derlc(DD,E,h3,k);
		prprogress(fid,' %6.4f   %6.3e\n',h3,F3);
		F1 = F2; F2 = F3;
		h1 = h2; h2 = h3;
		alf = alf*0.99; % decrease step size
	end
	h = h2;
	prwaitbar(0);

return

function F = derlc(DD,E,h,k) % crisp version
	% computation of the likelihood derivative for Parzen density
	% given distances D and their object minima E (for increased accuracy)
	m = size(DD,1);
	warning off MATLAB:divideByZero;
		Y = (DD-repmat(E,m,1))/(2*h*h); % correct for minimum distance to save accuracy
	warning on MATLAB:divideByZero;
	IY = find(Y<20);                % take small distance only, others don't contribute
	P = zeros(m,m);
	P(IY) = exp(-Y(IY));
	PP = sum(P,2)';
	FU = repmat(realmax,1,m);
	J = find(PP~=0); 
	FU(J) = 1./PP(J);
	FF = sum(DD.*P,2);
	warning off MATLAB:divideByZero;
		F = (FU*FF)./(h*h) - m*k;
	warning on MATLAB:divideByZero;
return

function h = parzenmls(A,fid) %soft version

	SS = gettargets(setlabtype(A,'soft'));
	[m,k,c] = getsize(A);
	DD= distm(+A) + diag(1e70*ones(1,m));
	E = min(DD);
	h = zeros(c,1);
	h0 = sqrt(max(E));    % initial estimate of h
	
	
	s = sprintf('parzenml: runover classes');
	prwaitbar(c,s,m > 100);
	iter = 0;
	
	for j=1:c
		prwaitbar(c,j)
		S = SS(:,j);
		h1 = h0;
		F1 = derls(DD,E,h1,k,S); % derivative

	  prprogress(fid,'parzenml: class %i : \n',j);
		prprogress(fid,' %6.4f   %6.3e\n',h1,F1);
		if abs(F1) < 1e-70 
			h(j) = h1;
			prwarning(4,'jump out\n');
			break;
		end
	
		a1 = (F1+m*k)*h1*h1;
		h2 = sqrt(a1/(m*k));  % second guess
		F2 = derls(DD,E,h2,k,S); % derivative

		prprogress(fid,' %6.4f   %6.3e\n',h2,F2);
		if (abs(F2) < 1e-70) | (abs(1e0-h1/h2) < 1e-6) 
			h(j) = h2;
			prwarning(4,'jump out\n');
			break;
		end
	
		% find zero-point of derivative to optimize h^2
		% stop if improvement is small, or h does not change significantly
	
		
		prwaitbar(100,'parzenml: Optimizing smoothing parameter',m > 100)
		iter = 0;
		alf = 1;
		while abs(1e0-F2/F1) > 1e-4 & abs(1e0-h2/h1) > 1e-3 & abs(F2) > 1e-70
			iter = iter+1;
			prwaitbar(100,100-100*exp(-iter/10));
			h3 = (h1*h1*h2*h2)*(F2-F1)/(F2*h2*h2-F1*h1*h1);
			if h3 < 0 % this should not happen
				h3 = sqrt((F2+m*k)*h2*h2/(m*k));
			else
				h3 = sqrt(h3);
			end
			h3 = h2 +alf*(h3-h2);
			F3 = derls(DD,E,h3,k,S);
			prprogress(fid,' %6.4f   %6.3e\n',h3,F3);
			F1 = F2; F2 = F3;
			h1 = h2; h2 = h3;
			alf = alf*0.99; % decrease step size
		end
		prwaitbar(0)
		h(j) = h2;
	end
  prwaitbar(0)
return

function F = derls(DD,E,h,k,S) %soft version
	% computation of the likelihood derivative for Parzen density
	% given distances D and their object minima E (for increased accuracy)
	% S are the object weigths
	c = size(S,2);                  % number of classes
	m = size(DD,1);
	Y = (DD-repmat(E,m,1))/(2*h*h); % correct for minimum distance to save accuracy
	IY = find(Y<20);                % take small distance only, others don't contribute
	F = 0;
	for j=1:c
		P = zeros(m,m);
		P(IY) = exp(-Y(IY));
		PP = S(:,j)'*P';
		FU = repmat(realmax,1,m);
		J = find(PP~=0);  
		FU(J) = S(J,j)'./PP(J);
		K = find(S(:,j)==0);
		FU(K) = zeros(1,length(K));
		FF = (DD.*P)*S(:,j);
		F = F + (FU*FF)./(h*h);
	end
	F = F - sum(S(:))*k;
return