lans_ppsproj.m

模式识别工具包

%	lans_ppsproj	- Project test data onto PPS in data/latent space
%	
%	[py, px, <,result>]	= lans_ppsproj(pps,y,<options>)
%
%	_____OUTPUTS____________________________________________________________
%	py	projected dataspace co-ordinate 	D x N	(col vectors)
%		NOT valid for -proj prob
%	px	projected LATENT co-ordinate 		Q x N	(col vectors)
%	result	optional results depending on method		(cell)
%		result{1}	? x N
%
%	_____INPUTS_____________________________________________________________
%	pps	Probabilistic Principal Surface			(struct)
%	y	test data in data space				(col vectors)
%	options	Projection specific options			(string)
%		-proj		Projection method
%			prob	mean of posterior p(x|y)	Q x N
%				result{1} : mode		Q x N
%				result{2} : like		scalar
%				result{3} : max(R)		1 x N
%
%			tri	nearest triangular projection
%			grid	nearest latent grid projection
%			nn*	nearest latent node
%				result{1} : squared error	1 x N
%			* default
%
%		-px		extra processing on px		(string)
%			none*	Default
%			unit	normalized to unit length
%			prop	normalized to (2-|x|)*xunit not exceeding
%				PXTHRESH = 1.25
%
%	_____NOTES______________________________________________________________
%	- for demo, call function without parameters
%		- demo relies on lans_pps demo
%	- py only valid for Euclidean based measures
%	- '-prob' slightly better than '-nn'
%
%	_____SEE ALSO___________________________________________________________
%	lans_proj2sqr lans_proj2tri lans_ppspost
%
%	(C) 2000.03.03 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/

%	_____TO DO______________________________________________________________

function	[py,px,result]	= lans_ppsproj(pps,y,options)

PXTHRESH	= 1.25;

if nargin>0
%__________ REGULAR ____________________________________________________________
if nargin<3
	options=[];
end

[D,N]	= size(y);
proj	= paraget('-proj',options,'nn');

px_o	= paraget('-px',options,'none');

lx	= lans_md2lin(pps.x);
tmpidx	= ones(1,N);

if length(pps.x)>1
	hascenter	= (size(pps.x{2},2)==3);
else
	hascenter	= 0;
end

%_____	if sphere
if strcmp(pps.manifold,'3sphere')
	issphere	= 1;
	wraparound	= 1;
	options		= paraset('-clos',1,options);
else
	issphere	= 0;
	wraparound	= 0;
end

%_____	initialize
f      		= pps.f;
result{1}	= realmax*ones(1,N);


switch proj
	case 'prob'
		[R,like]= lans_ppspost(pps,y,options);	% R is M x N
		px	= lx*R;				% Q x N
		py	= [];

		%__________	mode
		if nargout==3
			[result{3},tmpidx]= max(R); 		% 1 x N
			result{1}	= lx(:,tmpidx);
			result{2}	= like;
		end
	case 'nn'
		% ignore pseudo center node
		lf	= lans_md2lin(pps.f);
		d2	= lans_dist(lf(:,1:end-hascenter),y,'-Euclidean2');	% M x N
		[mind2,tmpidx]	= min(d2);
		px		= lx(:,tmpidx);
		py		= lf(:,tmpidx);
		%__________	
		if nargout==3
			result{1}	= mind2;
		end
	case 'grid'
		if pps.Q==1
			lf	= lans_md2lin(pps.f);
			lx	= lans_md2lin(pps.x);
			[py,px,result]	= lans_proj2line(y,lf,lx,options);
			return;
		end
		
		%_____	all vertical grids
		for grid_x=1:pps.xdim{1}(1)
			if issphere
				[tpy,tpx,tresult]	= lans_proj2line(y,[f{2}(:,1) squeeze(f{1}(:,grid_x,:)) f{2}(:,2)],[pps.x{2}(:,1) squeeze(pps.x{1}(:,grid_x,:)) pps.x{2}(:,2)]);
				
			else
				[tpy,tpx,tresult]	= lans_proj2line(y,squeeze(f{1}(:,grid_x,:)),squeeze(pps.x{1}(:,grid_x,:)));
			end
			wgood		= find(tresult{1}<result{1});
			result{1}(wgood)= tresult{1}(wgood);
			py(:,wgood)	= tpy(:,wgood);
			px(:,wgood)	= tpx(:,wgood);
		end
		%_____	all horizontal grids
		for grid_y=1:pps.xdim{1}(2)
			[tpy,tpx,tresult]	= lans_proj2line(y,squeeze(f{1}(:,:,grid_y)),squeeze(pps.x{1}(:,:,grid_y)),options);
			wgood		= find(tresult{1}<result{1});
			result{1}(wgood)= tresult{1}(wgood);
			py(:,wgood)	= tpy(:,wgood);
			px(:,wgood)	= tpx(:,wgood);
		end
	case 'tri'
		if pps.Q==1
			[py,px,result]	= lans_ppsproj(pps,y,paraset('-proj','grid',options));
			return;
		end
		hnum	= size(pps.x{1},2);
		vnum	= size(pps.x{1},3);
		%_____	grid patches
		for grid_y=1:(vnum-1)
			for grid_x=1:(hnum-1+wraparound)
				nextx		= mod(grid_x+1,hnum+1);
				nextx		= nextx+(nextx==0);
				fsqrpatch	= [f{1}(:,[grid_x nextx],grid_y) f{1}(:,[grid_x nextx],grid_y+1)]; 
				xsqrpatch	= [pps.x{1}(:,[grid_x nextx],grid_y) pps.x{1}(:,[grid_x nextx],grid_y+1)]; 
				[tresult,tpy,tpx]=lans_proj2sqr(y,fsqrpatch,xsqrpatch,'triangle');
				wgood		= find(tresult<result{1});
				result{1}(wgood)= tresult(wgood);
				py(:,wgood)	= tpy(:,wgood);
				px(:,wgood)	= tpx(:,wgood);
			end
		end
		
		%_____	non-grid patches (if applicable)
		if issphere
			p	= 0;
			for pole=[1 vnum];
				p=p+1;
				for grid_x=1:hnum
					nextx		= mod(grid_x+1,hnum+1);
					nextx		= nextx+(nextx==0);
					ftripatch	= [f{2}(:,p) f{1}(:,[grid_x nextx],pole)];
					xtripatch	= [pps.x{2}(:,p) pps.x{1}(:,[grid_x nextx],pole)];
					[tresult,tpy,tpx]	= lans_proj2tri(y,ftripatch,xtripatch);	
					wgood		= find(tresult<result{1});
					result{1}(wgood)= tresult(wgood);
					py(:,wgood)	= tpy(:,wgood);
					px(:,wgood)	= tpx(:,wgood);
				end	
			end
			
		end
	otherwise
		error('Unknown projection method.');
end
switch px_o
	case 'unit'
		px	= lans_unitvec(px);
	case 'prop'
		Q	= size(px,1);
		pxfac	= 2-sqrt(sum(px.*px));
		pxfac	= min(pxfac,PXTHRESH);
		pxfac	= ones(Q,1)*pxfac;
		px	= pxfac.*lans_unitvec(px);
end
%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_ppsproj.m in demo mode');

data_o  = '-linear 1 -repeat 0';
pps_o   = '-manifold 3sphere -alpha 1 -debug 0 -L1 5 -Lfac2 -Lratio 0 -regularize .00 -eint 30 -rint 60';
eint	= paraget('-eint',pps_o);
rint	= paraget('-rint',pps_o);
y       = lans_sphere([],[],[data_o ' ' pps_o]);

pps     = lans_ppsinit(y,pps_o);
y	= y + .5*randn(size(y));

%__________	plot manifold and data
clf;
lans_ppsplot(pps,'go-',[],y,'k.');

%__________	NN proj
[py,px,result]	= lans_ppsproj(pps,y,'-proj nn');
h(1)		= lans_plotseg(py,y,'r-','-hold 1');
leg{1}		= sprintf('MSE_{NN} = %1.3f',mean(result{1}));
nnresult	= result{1};

%__________	grid proj
[py,px,result]	= lans_ppsproj(pps,y,'-proj grid');
h(2)		= lans_plotseg(py,y,'b-','-hold 1');
leg{2}		= sprintf('MSE_{grid} = %1.3f',mean(result{1}));
gridresult	= result{1};


%__________	tri proj
[py,px,result]	= lans_ppsproj(pps,y,'-proj tri');
h(3)		= lans_plotseg(py,y,'k-','-hold 1');
leg{3}		= ['MSE_\Delta = ' sprintf('%1.3f',mean(result{1}))];
triresult	= result{1};

%__________	plot px (valid since latent and data space are equal)
h(4)		= lans_plotmd(px,'ko','-hold 1');
leg{4}		= 'px_\Delta';

%__________	prob proj
%[py,px,result]	= lans_ppsproj(pps,y,'-proj prob');
%lans_plotseg(px,y,'k-','-hold 1');


legend(h,char(leg),1);
%__________ DEMO ends __________________________________________________________
end