bsolk.m

模式识别工具包

%	bsolk		- Batch self-organizing line (k epochs)
%
%	function [nk,nf,trmse,vamse,temse] = bsolk(trdata,vadata,tedata,para,spara,k,f)
%
%	INPUTS
%	======
%	trdata	: training data 				(col vectors)
%	vadata	: validation data				(col vectors)
%		  =[] if none 
%	tedata	: test data					(col vectors)
%		  =[] if none
%	para	: see getpara.m					(row vector)
%	spara	: see getspara.m				(row vector)
%
%	*k	: starting principal curve knots		(row vector)
%	*f	: starting principal curve			(col vectors) 
%
%	OUTPUTS
%	=======
%	nk	: new principal curve knots			(row vector)
%	nf	: new principal curve				(col vectors)
%	trmse	: mse for training data				(scalar)
%	vamse	: mse for validation data 			(scalar)
%	temse	: mse for test data				(scalar) 
%
%	(C) 2000.06.28 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/

function [nk,nf,trmse,vamse,temse] = pcurvek(trdata,vadata,tedata,para,spara,k,f)

getpara;
getspara;

%----------	determine span value
spsize	= spdec*span;
spminv	= spmin*span;

%----------	create initial curve if none specified
if nargin<6,		%no initial principal curve specified
	[k,f]	= pinit(trdata,init,clos);
end

[kd,kn]	= size(k);
[fd,dum]= size(f);

%----------	create window buffer
win_trmse	= rvmake(win,1,1);
win_vamse	= rvmake(win,1,1);
win_temse	= rvmake(win,1,1);
win_vatol	= rvmake(win,1,1);
win_span	= rvmake(win,1,1);
win_k		= rvmake(win,kd,kn);
win_f		= rvmake(win,fd,kn);

oldvamse	= 1;

%----------	Compute initial errors
trmse	= projdist(trdata,k,f,para);

haste	= 0;
hasva	= 0;
if ~isempty(vadata),
	vamse 	= projdist(vadata,k,f,para);
	hasva	= 1;
else
	vamse	= trmse;
end

if ~isempty(tedata), 
	temse	= projdist(tedata,k,f,para);
	haste	= 1;
else
	temse	= trmse;
end

vatol	= 0;

%----------	Display Initial Curve
ts	= sprintf('[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f',0,trmse,vamse,temse,vatol,span);
if showg>0
	disp(ts);
end
if showg>1
	scatline(trdata,f,clos,ts);
	if showg>2
		axis equal;
	end
	drawnow;
end

%----------	do not exit within first win# iterations
for i	= 1:win
	[k,f,trmse]	= bsol1(trdata,para,k,f);
	if hasva==1
		vamse	= projdist(vadata,k,f,para);
	else
		vamse	= trmse;
	end
	if haste==1
		temse	= projdist(tedata,k,f,para);
	else
		temse	= trmse;
	end

	vatol		= abs(vamse-oldvamse)/oldvamse;
	oldvamse	= vamse;

	%----------	reduce span if validation error tolerance < sptol
	if vatol<sptol
		span	= max(spminv,span - spsize);
		putpara;
	end

	%----------	collect data over a window of size win
	win_trmse	= rvshift(trmse,win_trmse);
	win_vamse	= rvshift(vamse,win_vamse);
	win_temse	= rvshift(temse,win_temse);
	win_vatol	= rvshift(vatol,win_vatol);
	win_span	= rvshift(span,win_span);

	win_k		= rvshift(k,win_k);
	win_f		= rvshift(f,win_f);

	ts	= sprintf('[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f',i,trmse,vamse,temse,vatol,span);
	if showg>0
		disp(ts);
	end
	if showg>1
		scatline(trdata,f,clos,ts);
		if showg>2
			axis equal;
		end
		drawnow;
	end
end

broken	= 0;

%----------	after first win# iterations (window formed)
for i	= (win+1):iter
	[k,f,trmse]	= bsol1(trdata,para,k,f);

	if hasva==1
		vamse	 	= projdist(vadata,k,f,para);
	else
		vamse	= trmse;
	end
	if haste==1
		temse	= projdist(tedata,k,f,para);
	else
		temse	= trmse;
	end

	vatol		= abs(vamse-oldvamse)/oldvamse;
	oldvamse	= vamse;

	%----------	reduce span if validation error tolerance < sptol
	if vatol<sptol
		span	= max(spminv,span - spsize);
		putpara
	end

	%----------	collect data over a window of size win
	win_trmse	= rvshift(trmse,win_trmse);
	win_vamse	= rvshift(vamse,win_vamse);
	win_temse	= rvshift(temse,win_temse);
	win_vatol	= rvshift(vatol,win_vatol);

	win_k		= rvshift(k,win_k);
	win_f		= rvshift(f,win_f);

	ts	= sprintf('[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f',i,trmse,vamse,temse,vatol,span);
	if showg>0
		disp(ts);
	end
	if showg>1
		scatline(trdata,f,clos);
		title(ts)
		if showg>2
			axis equal;
		end
		drawnow;
	end

	%----------	Terminate if
	%			(a) vamse trend increases over a window win
	%			or
	%			(b) change in vamse dips below mtol

	trend	= win_vamse(2:win,2)-win_vamse(1:win-1,2);
	if sum(trend>=0)==(win-1)
		broken	= 1;
		break;
	elseif vatol<mtol
		break;
	end
end


if broken
	% the first record stored in window is the best before CV error inc.
	nk	= rvget(1,win_k,kd);
	nf	= rvget(1,win_f,fd);
	trmse	= rvget(1,win_trmse,1);
	vamse	= rvget(1,win_vamse,1);
	temse	= rvget(1,win_temse,1);
	vatol	= rvget(1,win_vatol,1);
	span	= rvget(1,win_span,1);
	siter	= i-win;
else
	% the iteration loop was exited normally, use last iterated value
	siter	= i;
	nk	= k;
	nf	= f;
end

ts	= sprintf('\n->[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f\n',siter,trmse,vamse,temse,vatol,span);
if showg>0
	disp(ts);
end
if showg>1
	scatline(trdata,f,clos);
	title(ts)
	if showg>2
		axis equal;
	end
	drawnow;
end