lans_ppsinit.m

模式识别工具包

		L		= L1+1;
		p_cen		= lans_scale(1:L1,GRIDMIN,GRIDMAX);
		p_sig		= (p_cen(2)-p_cen(1))*ones(1,L1);
		if gtm==1
			p_cen	= p_cen*L1/(L1-1);	% yield GTM scaling
			p_var	= p_fac*(p_sig.^2);	% operates on variance
		else
			p_var = (p_fac*p_sig).^2;	% " on sigma
		end

	case '2square'
		%__________	latent vectors x (arranged on 2-D grid)
		x1	= lans_scale(1:M1,GRIDMIN,GRIDMAX);
		x2	= lans_scale(1:M2,GRIDMIN,GRIDMAX);
		[pps.x{1}(2,:,:) pps.x{1}(1,:,:)]	= meshgrid(x1,x2);
		pps.xdim= lans_cellsize(pps.x,1);

		M	= M1*M2;

		%__________	latent basis (not arranged)
		l1	= lans_scale(1:L1,GRIDMIN,GRIDMAX);
		l2	= lans_scale(1:L2,GRIDMIN,GRIDMAX);
		[l(2,:,:) l(1,:,:)]	= meshgrid(l1,l2);
		p_cen	= lans_md2lin(l);
		h_dist	= l(1,2,1)-l(1,1,1);
		v_dist	= l(2,1,2)-l(2,1,1);
		L	= size(p_cen,2)+1;
		p_sig	= min(h_dist,v_dist)*ones(1,L-1);
		p_var	= (p_fac*p_sig).^2;	

	case '3sphere'
	% gridx gridy follows Cartesian convention
	% (not scanning convention, i.e. gridx ~= row, gridy ~=col)
		%__________	latent vectors x (arranged on 2-D grid)
		% single vector at N/S pole to test rotational invariance
		%
		% x{1}(:,gridx,gridy) gives latent col vector @ (gridx,gridy)

		% make elevation (elist) and rotation (rlist) indices
		eint	= paraget('-eint',options,10);
		rint	= paraget('-rint',options,10);

		elist	= -90:eint:90;
		ne	= length(elist);

		rlist	= 0:rint:350;
		nr	= length(rlist);
		pps.x	= lans_sphere(elist,rlist,options);
		
		% begin obsolete
		if center
			pps.x{2}	= [pps.x{2} zeros(Q,1)];
		end
		% end obsolete

		pps.xdim= lans_cellsize(pps.x,1);
		
		%__________	latent basis (not stored topologically)
		% - latent basis can follow the latent node or regular spacing
		% - only 1 basis needed at N and S pole
		% - L1 = max # nodes in each direction (x,y)
		
		% follows latent node's ratio of # x and y nodes
		if Lratio
			fac		= ne/nr;
			if ne>nr
				L_rlist	= lans_scale(1:(1+round(L1/fac)),0,360);
				L_elist	= lans_scale(1:L1,elist(0),elist(end));
			else
				L_rlist	= lans_scale(1:(L1+1),0,360);
				L_elist	= lans_scale(1:round(L1*fac),elist(0),elist(end));
			end
		else
			L_rlist	= lans_scale(1:(L1+1),rlist(1),360);
			L_elist	= lans_scale(1:L1,elist(1),elist(end));
		end

		L_rlist	= L_rlist(1:end-1);

		[p_cen,p_sig]	= lans_sphere(L_elist,L_rlist,'-repeat 0 -linear 1');
		p_var		= (p_fac*p_sig).^2;
		L		= length(p_sig)+1;	% reset correct L

	case '3hemisphere'		% center NOT FUNCTIONAL
	% gridx gridy follows Cartesian convention
	% (not scanning convention, i.e. gridx ~= row, gridy ~=col)
		%__________	latent vectors x (arranged on 2-D grid)
		% single vector at N/S pole to test rotational invariance
		%
		% x{1}(:,gridx,gridy) gives latent col vector @ (gridx,gridy)

		% make elevation (elist) and rotation (rlist) indices
		eint	= paraget('-eint',options,10);
		rint	= paraget('-rint',options,10);

		elist	= 0:eint:90;
		ne	= length(elist);

		rlist	= 0:rint:350;
		nr	= length(rlist);
		pps.x	= lans_sphere(elist,rlist);
			% ------- begin obsolete
			if center
				pps.x{2}	= [pps.x{2} zeros(Q,1)];
			end
			% ------- end obsolete
		pps.xdim= lans_cellsize(pps.x,1);
		
		%__________	latent basis (not stored topologically)
		% - latent basis can follow the latent node or regular spacing
		% - only 1 basis needed at N and S pole
		% - L1 = max # nodes in each direction (x,y)
		
		% follows latent node's ratio of # x and y nodes
		if Lratio
			fac		= ne/nr;
			if ne>nr
				L_rlist	= lans_scale(1:round(L1/fac),rlist(0),rlist(end));
				L_elist	= lans_scale(1:L1,elist(0),elist(end));
			else
				L_rlist	= lans_scale(1:L1,rlist(0),rlist(end));
				L_elist	= lans_scale(1:round(L1*fac),elist(0),elist(end));
			end
		else
			L_rlist	= lans_scale(1:L1,rlist(1),rlist(end));
			L_elist	= lans_scale(1:L1,elist(1),elist(end));
		end

		[p_cen,p_sig]	= lans_sphere(L_elist,L_rlist,'-repeat 0 -linear 1');
		p_var		= (p_fac*p_sig).^2;
		L		= length(p_sig)+1;	% reset correct L
		%------------------------------------- svar==?
	otherwise
		error('Manifold unsupported!');
end	% case manifold

%__________	compute fixed activation and bias
if (Lalpha~=1)&(strcmp(manifold,'3sphere'))
	% specify single unit vector emanating from origin 
	for l=1:L-1
		ovec(:,1,l)	= p_cen(:,l);
	end
	bas	= lans_basis('Gaussian-clamped',p_cen,'variance',p_var,'alpha',Lalpha,'omanifold',ovec);
else
	% isotropic latent basis
	bas	= lans_basis('Gaussian-spherical',p_cen,'variance',p_var,'unormalized',1);
end

lx	= lans_md2lin(pps.x);
M	= size(lx,2);

p	= lans_pdf(lx,bas,1);
po	= lans_pdf(zeros(Q,1),bas,1);

%__________	compute initial W
if Q==3
	% assume 3sphere
	mu	= mean(y')';
	sig	= std(y')';
	switch lower(init)
	case 'data'
		if center
			ry	= [y mu];
		else
			ry	= y;
		end
		if size(ry,2)~=size(p,2)
			error(sprintf('Initialization using data Failed! # data %d  ~= # nodes %d',N,M));
		end
	case 'gaussian'
		if center	% last point is mean
			ry	= [mu*ones(1,M-1)+sig*ones(1,M-1).*randn(D,M-1) mu];
		else
			ry	= mu*ones(1,M)+sig*ones(1,M).*randn(D,M);
		end
	case 'permute'
		if N<M
			error('# nodes > # data points!');
		end
		rorder	= randperm(N);
		ry	= y(:,rorder(1:M));
	case 'manifold'
		% sphere in principal subspace (3-D)
		[pc_data,pvar,paxis]    = lans_pca(y);	
		lx	= lans_md2lin(pps.x);
		scale	= sqrt(pvar(1:3))';
		ry	= (paxis(:,1:3).*(ones(D,1)*scale))*lx+mu*ones(1,M);
	end	% switch lower(init)
elseif Q==2
	% find 1st,2nd PC (form a grid) in data spc
	x1		= lans_scale(1:M1,GRIDMIN,GRIDMAX);
	x2		= lans_scale(1:M2,GRIDMIN,GRIDMAX);
	[pcay,svar]	= lans_pcgrid(x1,x2,y);
	% map latent grid to 2-D principal subspc in data spc
	ry		= lans_md2lin(pcay);
elseif Q==1
	% find 1st PC in data spc
	[pcay,svar]	= lans_pcgrid(pps.x{1},y);
	% map latent line to 1st PC in data spc
	ry		= pcay;
end
W	= ry/p;
%_____________________________________________________________
pps.L		= L;
pps.M		= M;
pps.p		= p;
pps.po		= po;
pps.p_cen	= p_cen;
pps.p_sig	= p_sig;
pps.W		= W;

%	reference vectors in original space
lf		= W*p;
pps.f		= lans_lin2md(lf,pps.xdim);

%__________	Compute Jacobian gradient at each reference vector x __________
%at each basis p_cen,	dp/dx	= (p_cen-x)*p/(sig^2)
if (alpha~=1)|debug
	pps.dpdx	= zeros(L,Q,M);				% L x Q x M
%	[lx,xdim]	= lans_md2lin(pps.x);
	sigma		= ones(M,1)*(pps.p_sig.*pps.p_sig);	% M x L
	for q=1:Q
		muplane	= ones(M,1)*pps.p_cen(q,:);		% M x L
		xplane	= lx(q,:)'*ones(1,L-1);			% M x L
		dplane	= (muplane-xplane)./(sigma.*sigma);	% M x L
		gplane	= dplane.*p(1:end-1,:)';		% M x L
		pps.dpdx(1:L-1,q,:)	= gplane';		% L-1 x Q x M
		pps.dpdx(L,q,:)		= 0;
	end
end

%__________	initial beta (median dist between nodes in data space)
fdist2	= lans_dist(lf,'-metric Euclidean2') + eye(M)*realmax;
% mean not used here because it is skewed by the zero distances of
% duplicating pole values
%nnd2	= mean(min(fdist2));
nnd2	= median(min(fdist2));

if Q<3
	% also consider PCA eigenvalues
	pps.beta= 1/max(lans_clip(nnd2,eps,realmax),svar(Q+1));
else
	% use median distance between nodes in data space
	pps.beta= 1/lans_clip(nnd2,eps,realmax);
end

%__________	prior for reference vectors
pps.prior	= ones(1,M)/M;		% constant in general

%__________	initial manifold gradient & sigma 
if alpha~=1
	dfdx		= W*reshape(pps.dpdx,L,M*Q);
	pps.dfdx	= reshape(dfdx,D,Q,M);	% NOT orthogonal but spans R^Q
end

%__________	likelihood
[pps.R,Lcomp,Lavg]	= lans_ppspost(pps,y,options);


%__________	projection
if nargout>3
	[dum,dum,result]	= lans_ppsproj(pps,y,options);	
	mse			= mean(result{1});
end

%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_ppsinit.m in demo mode');
N	= 100;          % # samples
nsig	= 0;		% noise s.d.
%nsig	= .05;		% noise s.d.

%pps_o	= '-alpha 1.1 -center 1 -init data -Lfac 2 -L1 7 -Lratio 0 -Lalpha 1 -M1 10 -manifold 3sphere -eint 30 -rint 30';
pps_o	= '-alpha 1.1 -center 1 -init manifold -Lfac 2 -L1 7 -Lratio 0 -Lalpha 1 -M1 10 -manifold 3hemisphere -eint 30 -rint 30';
%pps_o	= '-alpha 1.1 -Lfac 2 -L1 10 -Lratio 1 -Lalpha 1 -M1 10 -manifold 2square';
%pps_o	= '-alpha 1.1 -Lfac 2 -L1 10 -Lratio 1 -Lalpha 1 -M1 10 -manifold 1line';

y	= lans_sphere(0:30:90,0:30:330,'-repeat 0 -linear 1');
y	= y+nsig*randn(size(y));
%y	= rand(3,N);

%__________     Initialize PPS
[pps,Lcomp,Lavg]	= lans_ppsinit(y,pps_o)

%__________	plot results
clf;
h{1}	= lans_ppsplot(pps,'bo-',pps_o,y,'g.');
h{2}	= lans_plotmd(pps.p_cen(:,[end-1 1:end-2 end]),'rx-','-hold 1');	% latent basis

h{3}	= lans_plotmd(pps.W*pps.po,'r+','-hold 1');
%_________	plot gradients
lf	= lans_md2lin(pps.f);
lans_plotmd(zeros(3,1),'ko','-hold 1');		% center
lans_plotmd(lf(:,end),'ko','-hold 1');	% north pole
lans_plotgrad(lf(:,end),pps.dfdx(:,:,end),'b')

hh	= cat(2,h{:})
title('R^{Q} : Latent Space (Rotate with mouse)');
legend(hh,'pps.f','y','pps.p\_cen',1);
rotate3d on;

%mse	= sum(vdist2(y,yhat))
%__________ DEMO ends __________________________________________________________
end