lans_pcurvek.m

模式识别工具包

%	lans_pcurvek	- Principal Curve iteration
%
%	[nk,nf,mse] = lans_pcurvek(data,para,k,f)
%
%	_____OUTPUT_____________________________________________________________
%	nk	new principal curve knots			(row vector)
%	nf	new principal curve				(col vectors)
%	mse	mean squared error				(col vector)
%		mse(1)	training mse				(scalar)
%		mse(2)	validation mse		-1 if none	(scalar)
%		mse(3)	test mse		-1 if none	(scalar)
%
%	_____INPUT______________________________________________________________
%	data	must specify all 3 types			[structure]
%		.train	training data				(col vectors)
%		.valid	validation data				(col vectors)
%		.test	test data				(col vectors)
%		[] if empty
%	para	see lanspara.m paraget.m			(string)
%		-algo -clos -datafile -init -iter -kernel -mtol
%		-resultfile -rsort -showg -showt -smoopara -smoother
%		-span -spdec -spmin -sptol -win
%		[] if empty
%	k	starting principal curve knots			(row vector)
%		[] if empty
%	f	starting principal curve			(col vectors) 
%		[] if empty
%
%	_____SEE ALSO___________________________________________________________
%	pcurve1 paraget
%
%	(C) 1998.04.20 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,mse] = pcurvek(data,para,k,f)

%----------	argument check
trdata	= data.train;

vadata	= data.valid;
hasva	= ~isempty(vadata);

tedata	= data.test;
haste	= ~isempty(tedata);

%----------	argument check	
if isempty(k),
	clos	= paraget('-clos',para);
	init	= paraget('-init',para);
	[k,f]	= pinit(trdata,init,clos);	% creat my own starting curve
end

%----------	get parameters	
iter	= paraget('-iter',para);
mtol	= paraget('-mtol',para);
showg	= paraget('-showg',para);
showt	= paraget('-showt',para);
span	= paraget('-span',para);
spdec	= paraget('-spdec',para);
spmin	= paraget('-spmin',para);
sptol	= paraget('-sptol',para);
win	= paraget('-win',para);

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

%----------	create initial projection (with corresponding mse)
[ptrdata,pk,pf,trmse]	= project(trdata,k,f);

[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
if hasva
	vamse 	= projdist(vadata,k,f,para);
else
	vamse	= trmse;
end

if haste 
	temse	= projdist(tedata,k,f,para);
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 showt==1
	disp(ts);
end

if showg==1
	scatline(trdata,f,0,para);
end

%========== do not exit within first win# iterations
for i	= 1:win
	[nk1,nf1,ptrdata,pk,pf,trmse]	= pcurve1(ptrdata,pk,pf,para);
	if hasva==1
		vamse	= projdist(vadata,nk1,nf1,para);
	else
		vamse	= trmse;
	end

	if haste==1
		temse	= projdist(tedata,nk1,nf1,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);
		para	= paraset('-span',span,para);
	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(nk1,win_k);
	win_f		= rvshift(nf1,win_f);

	ts	= sprintf('[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f',i,trmse,vamse,temse,vatol,span);
	if showt==1
		disp(ts);
	end
	if showg==1
		scatline(trdata,nf1,0,para);
	end
end

broken	= 0;

%----------	after first win# iterations (window formed)
%==========================================================
for i	= (win+1):iter
	[nk1,nf1,ptrdata,pk,pf,trmse]	= pcurve1(ptrdata,pk,pf,para);
	if hasva==1
		vamse	= projdist(vadata,nk1,nf1,para);
	else
		vamse	= trmse;
	end
	if haste==1
		temse	= projdist(tedata,nk1,nf1,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);
		para	= paraset('-span',span,para);
	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(nk1,win_k);
	win_f		= rvshift(nf1,win_f);

	ts	= sprintf('[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f',i,trmse,vamse,temse,vatol,span);
	if showt==1
		disp(ts);
	end
	if showg==1
		scatline(trdata,nf1,0,para);
	end

	%----------	Terminate if
	% (a) vamse trend increases over a window win
	%	or
	% (b) vamse(=trmse if no validation set) dips below mtol
	%

	trend	= win_vamse(2:win,2)-win_vamse(1:win-1,2);
	if (sum(trend>=0)==(win-1))&hasva
		broken	= 1;
		break;
	elseif vatol<mtol		% vatol = trtol if no validation set
		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);
	mse(1)		= rvget(1,win_trmse,1);
	mse(2)		= rvget(1,win_vamse,1);
	mse(3)		= 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		= nk1;
	nf		= nf1;
	mse(1)		= trmse;
	mse(2)		= vamse;
	mse(3)		= temse;
end

ts	= sprintf('\n->[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f\n',siter,trmse,vamse,temse,vatol,span);

str	= sprintf('[%3d]tr=%.4f va=%.4f te=%.4f tol=%.4f span=%0.2f',siter,trmse,vamse,temse,vatol,span);

if showt==1
	disp(ts);
end;

if showg==1
	scatline(trdata,nf,0,para);
end;