lans_project.m

模式识别工具包

%	lans_project	- Project data onto a principal curve
%
%	[prsurf,prdata,mse,se]	= lans_project(data,psurf[,para])
%
%	_____OUTPUTS____________________________________________________________
%	prsurf	projection of points onto current psurf.f	(structure)
%		see lans_psurf.m
%	prdata	ordered data w.r.t. prsurf.x			(col vectors)
%	mse	mean squared error of data to prsurf.pf		(scalar)
%	se	squared error/dist of each data to prsurf.pf	(row vector)
%
%	_____INPUTS_____________________________________________________________
%	data	data points (# may differ from #knots)		(col vectors)
%	psurf	principal curve					(structure)
%		see lans_psurf.m
%	para	see lanspara.m paraget.m			(string)
%		-algo -clos -rsort -smoother
%
%	_____NOTES______________________________________________________________
%	prec not much help in maintaining same output for MATCOM and MATLAB
%	errors accumulate
%
%	- mse, se are SQUARED errors, NOT distances!
%	- prsurf is a sorted and projected version of psurf
%	- psurf.x and psurf.f MUST be sorted!!!
%	- OK if # data points ~= # knots	
%
%	- if latent dimensions use different units specify psurf.xQ
%	  to get interpolated values in prsurf.xQ
%
%	_____SEE ALSO___________________________________________________________
%	lans_expect
%
%	(C) 1999.11.11 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	[prsurf,prdata,mse,se]	= lans_project(data,psurf,para,debug)

if nargin<3,
	para	= [];
	if nargin<1
		clc;
		help lans_project;
		return;
	end
end

%----- get parameters
algo	= paraget('-algo',para);
clos	= paraget('-clos',para);
smoother= paraget('-smoother',para);
rsort	= paraget('-rsort',para);

%-----	preassign prsurf=psurf to keep probability structures in prsurf
% this step is vital in the iteration of the EM based algorithm in lans_expect.m
if algo>2
	prsurf	= psurf;
end

%-----	cubic spline smoother needs sorted knots
%if smoother==3
%	rsort	= 1;
%end

[d N]	= size(data);		% N	= # of data points
[d M]	= size(psurf.f);	% M	= # of knots


%The following is used instead of ZERO due to difference in precision
%across platforms
%ALSO, in determining smallest distance, if difference between min knot
%and min segment to point is less than prec, we use the seg dist

%prec	= 1e-7;
prec	= 0;

k	= psurf.x;
f	= psurf.f;
%---------- compute length of closing segment
if clos~=0
	clos		= vdist(f(:,M),f(:,1));
	f(:,M+1)	= f(:,1);
	k(M+1)		= k(M)+clos;
	M2		= M + 1;
else
	M2		= M;
end

bound	= zeros(3,N);	% lower/upper bounds and alpha for each projected point

%---------- project points onto nearest interval to form M,nf
for p	= 1:N
	x	= data(:,p);			% current point p
	rx	= x*ones(1,M2);			% replicate x M times

	%sdist2	= 1e5;

	%-----	Compute dist to each knot (in order)
	k2x	= rx-f;				% vector: all knots to x 1:M
	kdist2	= vdist2(k2x);			% sqr dist to all knots	 1:M		
	
	%-----	Find closest KNOT!!!
	minkdist2= min(kdist2(1:M));		% smallest dist to knot

	%-----	Compute unit segment vector
	l2r	= f(:,2:M2)-f(:,1:M2-1);	% compute segment vector
	dl2r	= vdist(l2r);			% compute segment length
	if dl2r>0
		l2r1	= l2r./(ones(d,1)*dl2r);% unit segment vector
	else
		l2r1	= l2r;
	end

	%-----	Compute projections onto each segment and find valid ones
	projd	= sum(k2x(:,1:M2-1).*l2r1);	% projected dist on all segments

	vidx1	= find(projd>prec);			% keep ones on right
	vidx2	= find(projd(vidx1)<dl2r(vidx1));	% keep those on segment
	vidx	= vidx1(vidx2);

	if ~isempty(vidx)
		% nearest lies on a segment
		onseg	= (ones(d,1)*projd(vidx)).*l2r1(:,vidx); % proj. vector
		onsegv	= onseg+f(:,vidx); % pos of proj. vector w.r.t. origin
		dist22x	= vdist2(onsegv,x*ones(1,length(vidx)));
		minsdist2= min(dist22x);
		minsidx	= min(find(dist22x==minsdist2));
		if minsdist2<minkdist2
			% on segment
			pf(:,p)	= onsegv(:,minsidx);
			idx	= vidx(minsidx);
			pk(p)	= k(idx)+projd(idx);
			alpha	= projd(idx)/dl2r(idx);
			bound(:,p)=[idx;idx+1;alpha];
		else
			% knot is nearer
			minkidx	= min(find(kdist2==minkdist2));
			pf(:,p)	= f(:,minkidx);
			pk(p)	= k(minkidx);
			bound(:,p)=[minkidx;minkidx;0]; 
		end
	else
		% nearest is the knot nearestknot
		minkidx	= min(find(kdist2==minkdist2));
		pf(:,p)	= f(:,minkidx);
		pk(p)	= k(minkidx);

		bound(:,p)	= [minkidx;minkidx;0];
	end
	
end	%p

%---------- sort the new k/f
%pk	= sqrt(pk);

prsurf.x	= pk;
prsurf.f	= pf;

if isfield(psurf,'xQ')
	Q	= size(psurf.xQ,1);
	lidx	= bound(1,:);
	ridx	= bound(2,:);
	a	= ones(Q,1)*bound(3,:);
	a1	= 1-a;
	prsurf.xQ=psurf.xQ(:,lidx).*a1 + psurf.xQ(:,ridx).*a;
end

if rsort
	[pk,ind]= sort(pk);
	%-----	sort corresponding probabilistic structures
	prsurf	= lans_porder(prsurf,ind,para);
	prdata	= data(:,ind);
else				% no reorder
	prdata	= data;
end

se	= vdist2(prdata,prsurf.f);
mse	= mean(se')';