lans_psurfinit.m

模式识别工具包

%	lans_psurfinit	- Principal surface initialization/creation (GTM based)
%
%	[psurf,prsurf,mse,se]	= lans_psurfinit(y)
%				= lans_psurfinit(y,options)
%
%	_____OUTPUTS____________________________________________________________
%	psurf	principal surface structure			(structure)
%		.f	psurf values f(x) in original dimension	(col vectors)
%		.D	dimensionality of original space	(integer+)
%		.covars	covariance of y about .f		(2/3-D array)
%			for GTM, .covars contain current
%			covars (which may be oriented)
%		.covartype	covariance type			(string)
%				'spherical'	.covars = row vector	(1xM)
%				'diagonal'	.covars = col vectors	(DxM) 
%				'full'		.covars = 3-D array	(DxDxM)
%		.priors	prior probability of each knot		(row vector)
%
%		----- GTM -----
%		.beta	GTM inverse spherical variance		(scalar)
%		.lact	GTM latent basis activation w.r.t. .x	( (L+1)xM )
%			row j represents the activation at
%			latent reference vector j
%		.lbas	GTM latent basis centers (QxL)		(col vectors)
%		.lgrad	GTM latent basis gradient w.r.t. x	( LxQxM )
%			row (:,:,j) represents the latent basis
%			gradient at reference vector j where the
%			LxQ matrix is the Jacobian
%			dphi/dx	= [
%				dphi_1(x)/dx_1 ... dphi_1(x)/dx_Q
%
%				dphi_L(x)/dx_1 ... dphi_L(x)/dx_Q
%				= [dPhi/dx_1 ... dPhi/dx_Q]
%
%			dphi_l(X)/dX=(Mu_l-X)phi_l(X)/(s^2)
%					
%			Gradient vectors in original dimension can be obtained
%			via
%			df/dq	=psurf.lgrad(:,q,:)*psurf.W(1:end-1,:);
%			which is a MxDxQ matrix
%			* lgrad is only approximately orthogonal!!!!
%			* resulting in an approximately orthogonal grad!!!
%		.L	GTM # of latent basis functions		(integer+)
%			must be k^2 for 2square manifold
%		.M	# of knots/latent vectors		(integer+)
%			must be k^2 for 2square manifold
%		.Q	GTM dimensionality of manifold		(integer+)
%		.s	GTM latent basis width 			(scalar)
%			(fraction to nearest basis)
%			for latent basis
%		.W	GTM Weight matrix			( Dx(L+1) )
%		.x	GTM knots in latent/lower dimension	(col vectors)
%		.xidx	manifold indices to .x and .f		(Q-D vector)
%			e.g. .xidx(3,4) will denote index of ref vector
%			@ latent manifold row 3, col 4
%		Internal fields
%		---------------
%		.Cinv	inverses of the oriented .covars	(DxDxM)
%		.Cdet	squared determinant of oriented .covars	(1xM)
%		.Cun	unoriented covars			(DxDxM)
%	prsurf	subset of psurf					(psurf)
%		see lans_psurfproj.m
%	mse	Mean squared error of y to prsurf		(scalar)
%		USES Distance to nearest node aproximation
%	se	squared error/dist of each data to psurf.f	(row vector)
%		USES Distance to nearest node aproximation
%
%	_____INPUTS_____________________________________________________________
%	y	unormalized data in original dimension		(col vectors)
%	options	see lanspara.m					(string)
%		GTM specific options
%		beta	inverse variance			(scalar)
%		manifold					(string)
%		L	# of latent bases			(integer)
%		s	standard deviation factor w.r.t.	(scalar)
%			nearest latent basis e.g.
%			1 => s.d. = 1xwidth to nearest neighboring basis
%		
%
%	_____NOTES______________________________________________________________
%	- Use lans_psurfproj to get better projection estimate
%	- Needs GTM toolbox version 1.02 or higher
%	- covars are oriented here if specified
%	- uses randn.m, seed randn('seed',n);
%	- set beta = 0 to use algorithm initialized values, otherwise
%	specify beta. If unspecified, a default value of 1 will be used.
%	- code added to ensure that first dimension of eigenvector is >0 (+ve)
%	  this is to take care of the eigs.m discrepencies between matlab and
%	  matcom
%	- * Note that each psurf.lgrad is only approximately orthogonal!!!!
%
%	_____SEE ALSO___________________________________________________________
%	lans_psurfproj lans_psurfgrad
%
%	(C) 1999.07.13 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 [psurf,prsurf,mse,se]= lans_psurfinit(y,options)
if nargin<2
	options	= [];
	if nargin<1
		clc;
		help lans_psurfinit;
		break;
	end
end

algo		= paraget('-algo',options);
beta		= paraget('-beta',options);	% initial inverse variance
clos		= paraget('-clos',options);
manifold	= paraget('-manifold',options);
orient		= paraget('-orient',options);
psurf.covartype	= paraget('-covartype',options);
pinittype	= paraget('-pinittype',options);

[D N]		= size(y);
psurf.D		= D;
psurf.L		= paraget('-L',options);
psurf.M		= paraget('-M',options);
	if psurf.M==0
		psurf.M	= N;
	end
psurf.s		= paraget('-s',options);

%__________	Initialize manifold
switch (manifold)
	case '1line'
		psurf.Q	= 1;
		if psurf.L==0
			psurf.L	= ceil(psurf.M/psurf.s);
		end
		[xT,lbasT,lactT,WT,psurf.beta]	= gtm_stp1(y',psurf.M,psurf.L,psurf.s);
		%	in col vector notation
		psurf.lact	= lactT';
		psurf.lbas	= lbasT';
		psurf.W		= WT';
		psurf.x		= xT';
		%	sigma	= latent basis' width 
		sigma		= psurf.s*(psurf.lbas(2)-psurf.lbas(1));
	case '2square'		%***************** NEEDS WORK
		psurf.Q	= 2;
		%_____	make sure sqrt(psurf.M) is an interger
		Mup	= ceil(sqrt(psurf.M));
		Mlow	= floor(sqrt(psurf.M));
		if Mup~=Mlow
			if Mup<N
				psurf.M	= Mup*Mup;
			else
				psurf.M	= Mlow*Mlow;
			end
		end
		%_____	make sure sqrt(psurf.L) is an interger
		Lup	= ceil(sqrt(psurf.L));
		Llow	= floor(sqrt(psurf.L));
		if Lup~=Llow
			if Lup<sqrt(psurf.M)
				psurf.L	= Lup*Lup;
			else
				psurf.L	= Llow*Llow;
			end
		end
		[xT,lbasT,lactT,WT,psurf.beta]	= gtm_stp2(y',psurf.M,psurf.L,psurf.s);
		%	in col vector notation
		psurf.lact	= lactT';
		psurf.lbas	= lbasT';
		psurf.W		= WT';
		psurf.x		= xT';
		%	sigma	= latent basis' width 
		sigma		= psurf.s*(psurf.lbas(2,1)-psurf.lbas(2,2));
		psurf.xidx	= reshape(1:psurf.M,[sqrt(psurf.M) sqrt(psurf.M)]);

	case '3sphere'		%********************** NEEDS WORK
		psurf.Q	= 3;
		error('Unimplemented');
	otherwise
end

%	reference vectors in original space
psurf.f		= psurf.W*psurf.lact;

%__________	Compute Jacobian gradient at each reference vector x
for q=1:psurf.Q
	% row q of L latent basis centers repeated M times
	% to give MxL center plane
	muplane		= ones(psurf.M,1)*(psurf.lbas(q,:));
	% L cols of repeated ref vectors' q values 
	% row(latent dim) q of M samples repeated L times
	% to give MxL sample plane
	xplane		= psurf.x(q,:)'*ones(1,psurf.L);
	% the difference plane MxL
	dplane		= (muplane-xplane)/(sigma*sigma);
	% MxL grad plane for dimension q
	gradplane		= dplane.*lactT(:,1:end-1);
	psurf.lgrad(:,q,:)	= gradplane';
end
%__________	reset beta to user specified value (if not 0)
if beta~=0
	psurf.beta	= ones(1,psurf.M)*beta;
else
	beta		= ones(1,psurf.M)*psurf.beta;
end

%----------	set up mixture parameters (needed later for sorting)
psurf.priors	= ones(1,psurf.M)/psurf.M;
if algo==4	% GTM uses only spherical covariances		
	psurf.covars	= ones(1,psurf.M)./beta;
elseif algo==3
	switch psurf.covartype
	case 'spherical'
		psurf.covars	= ones(1, psurf.M)./beta;
	case 'diagonal'
    		psurf.covars	= ones(D, psurf.M)./beta;
    	case 'full'
   		psurf.covars	= repmat(eye(D), [1 1 psurf.M])./beta;  
	end
end

%__________	orient covariances and put them in .covars if specified
if orient&(algo>2)
	psurf	= lans_orient(psurf,options);
end

%----------	compute data projection onto psurf and corresponding MSE
[prsurf,mse,se]		= lans_psurfproj(y,psurf);