lans_psurf1.m

模式识别工具包

%	lans_psurf1	- 1 iteration of psurf (GTM) training
%
%	[psurf]	= lans_psurf1(y,psurf,options)
%
%	_____OUTPUTS____________________________________________________________
%	psurf	principal surface structure			(structure)
%		see lans_psurfinit.m
%	like	current likelihood				(scalar)
%		
%	_____INPUTS_____________________________________________________________
%	y	unormalized data in original dimension		(col vectors)
%	psurf	principal surface structure			(structure)
%		see lans_psurfinit.m
%	options	see lanspara.m					(string)
%		GTM specific options
%		beta	inverse variance			(scalar)
%		L	# of latent bases			(integer)
%		s	standard deviation factor w.r.t.	(scalar)
%			nearest latenet basis e.g.
%			1 => s.d. = 1xwidth to nearest neighboring basis
%		
%
%	_____NOTES______________________________________________________________
%	for demo, call function without parameters
%
%	- ASSUMES that covars is properly oriented (if desire orientation)
%	- This is a modified version of gtm_trn.m from the GTM toolbox
%	- Needs GTM toolbox version 1.02 or higher
%	- when orientation is used, psurf.covars is used to store the modified
%	  oriented covariances, and psurf.beta will store the single common
%	  isotropic variance value beta
%
%	D	original dimensionality
%	Q	reduced dimensionality (of manifold)
%	L	# of latent basis
%	L1	L + 1 (bias term)
%	M	# of reference vectors/nodes
%	N	# of samples
%
%	_____SEE ALSO___________________________________________________________
%	lans_psurfinit	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,like]	= lans_psurf1(y,psurf,options)

if nargin>0
%__________ REGULAR ____________________________________________________________
%mode		= paraget('-mode',options);
mode		= 0;
attenuate	= paraget('-attenuate',options);
regularize	= paraget('-regularize',options);
orient		= paraget('-orient',options);
if attenuate==1
	orient	= 0;
end

covarmax	= paraget('-covarmax',options);
covarmin	= paraget('-covarmin',options);
covarclip	= paraget('-covarclip',options);

WT		= psurf.W';
beta		= psurf.beta;
if length(beta)>1
	beta	= beta(1);
end
yT		= y';
[D N]		= size(y);
actT		= psurf.lact;
act		= actT';		%	M	= # latent ref. vectors
[M,L1]		= size(act);		%	L1	= # latent basis + 1
ND		= N*D;

% Declare global variables
global gtmGlobalDIST;
global gtmGlobalR;
global gtmGlobalRfactor;		% added by Kui-yu

if orient
	gtmGlobalRfactor	= psurf.Cidet2*ones(1,N);
end

if (mode > 0)
	global gtmGlobalMinDist;
	global gtmGlobalMaxDist;
end

% Pre-allocate matrices
gtmGlobalDIST	= zeros(M, N);
gtmGlobalR	= zeros(M, N);
if (mode > 0)
	gtmGlobalMinDist = zeros(1, N);
	gtmGlobalMaxDist = zeros(1, N);
end
A	= zeros(L1, L1);
cholDcmp= zeros(L1, L1);

% Use a sparse representation for the weight regularizing matrix.
if (regularize > 0)
	LAMBDA		= regularize*speye(L1);
	LAMBDA(L1, L1)	= 0;
end  

%-----------------------------------------------------------------
% Calculate initial distances (alters global : gtmGlobalDIST)
% in DATA space between reference vectors and data
%
% gtmGlobalDIST	is M x N
refrvecs	= act*WT;		% reference row vectors
if orient==0
	gtm_dstg(yT, refrvecs, mode);
else
	gtmGlobalDIST	= lans_mahadist(refrvecs',y,psurf.Cinv);
	% The above has not been weighted by the determinant
end

%-----------------------------------------------------------------
%Training
% Calculate responsabilities (alters global : gtmGlobalDIST, gtmGlobalR
like	= gtm_rspg1(beta, D, mode, orient);

% Calculate matrix be inverted (act'*G*act + lambda*I in the papers).
% Sparse representation of G normally executes faster and saves
% memory

%gtmGlobalR(:,68:71)	%//////////////////////////////////////////////////

if regularize>0
	A = full(actT*spdiags(sum(gtmGlobalR')', 0, M, M)*act + (LAMBDA./beta));
else
	A = full(actT*spdiags(sum(gtmGlobalR')', 0, M, M)*act);
end

% A is a symmetric matrix likely to be positive definite, so try
% fast Cholesky decomposition to calculate WT, otherwise use SVD.
% (actT*(gtmGlobalR*yT)) is computed rigth-to-left, as gtmGlobalR
% and yT are normally (much) larger than actT.
[cholDcmp singular] = chol(A);
if (singular)
	WT = pinv(A)*(actT*(gtmGlobalR*yT));
else
	WT = cholDcmp \ (cholDcmp' \ (actT*(gtmGlobalR*yT)));
end

% Assign new Weight Matrix
psurf.W	= WT';

% Assign new reference vectors
refrvecs	= act*WT;		% reference row vectors
psurf.f		= refrvecs';

%-----------------------------------------------------------------
% Calculate new distances using updated W (alters global : gtmGlobalDIST)
if orient==0
	gtm_dstg(yT, refrvecs, mode);
else
	gtm_dstg(yT, refrvecs, mode);
%	For oriented covars, the above is an approximated M-step
%	gtmGlobalDIST	= lans_mahadist(refrvecs',y,Cinv);
%	The above has not been weighted by the determinant
end

% Calculate new value for beta scalar
beta	= ND / sum(sum(gtmGlobalDIST.*gtmGlobalR));
if covarclip
	psurf.beta	= 1/lans_clip(1/beta,covarmin,covarmax);
else
	if (beta>covarmin)&(beta<covarmax)
		psurf.beta	= beta;
	end
end

psurf.beta	= ones(1,M)*beta;

% Assign covariances if not oriented
if orient==0
	psurf.covars	= 1./psurf.beta;
else
	psurf		= lans_orient(psurf,options);
end

clear global gtmGlobalDIST gtmGlobalR gtmGlobalRfactor gtmGlobalMinDist gtmGlobalMaxDist;


return

%_____	The following fails
rsum	= sum(gtmGlobalR');
for m=1:psurf.M
	d2	= vdist2(y,psurf.f(:,m)*ones(1,N))';
	svar(m)	= gtmGlobalR(m,:)*d2/(D*rsum(m));
	if covarclip
		beta(m)	= 1/lans_clip(svar(m),covarmin,covarmax);
	else
		if (svar(m)>covarmin)&(svar(m)<covarmax);
			beta(m)	= 1/svar(m);
		end
	end
end
%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_psurf1.m in demo mode');
%__________	Parameters
N		= 100;		% # samples
K		= 5;		% # iterations
dstring = '-algo 4 -attenuate 1 -beta 0 -L 2 -M 25 -manifold 2square -orient 1 -s 2 -tol 0.001 ';

%__________	Generate memispherical data
[y,ball]	= lans_genball(N);
ball.cons	= 'hemisphere';
[y,ball]	= lans_genball(N,0,ball);

%__________	Initialize psurf
[psurf,prsurf,mse,se]   = lans_psurfinit(y,dstring);

%__________	Iterate
for k=0:K
	if k>0
		psurf	= lans_psurf1(y,psurf,dstring);
		if k==1
			title('Initial Surface - Press any key to continue')
			pause
		end
	end
	clf;
	lans_psurfpgrid(psurf);	% draws grid
	grad	= lans_psurfgrad(psurf);
	lans_plotgrad(psurf.f,grad);
	title(k);
	drawnow;
end

%__________ DEMO ends __________________________________________________________
end