changeset 15:a6d5597bd922

Trying alternative fixed point computation in mc_global_info1>mc_fixpt2;
Also factored out randnat from some other functions.

--- a/mt_get_transmat_at.m	Sun Feb 26 23:17:29 2012 +0000
+++ b/mt_get_transmat_at.m	Mon Feb 27 14:41:24 2012 +0000
@@ -22,7 +22,7 @@
 	% distance from target
 	d2 = sum((info(:,1:2) - repmat(target,L,1)).^2,2);
 	M = find(d2==min(d2));
-	J = M(1+floor(length(M)*rand));
+	J = M(randnat(length(M)));
 	T = tmats(:,:,J);
 	I = info(J,:);
 	if Sys.shuffle

--- a/mt_get_transmat_near.m	Sun Feb 26 23:17:29 2012 +0000
+++ b/mt_get_transmat_near.m	Mon Feb 27 14:41:24 2012 +0000
@@ -29,7 +29,7 @@
 
 	M = find(Mask);
 %	J = M(argmax(info(M,3)));
-	J = M(1+floor(length(M)*rand));
+	J = M(randnat(length(M)));
 	T = tmats(:,:,J);
 	I = info(J,:);
 	%Mask=2*int16(Mask); Mask(J)=1;

--- a/mt_resample.m	Sun Feb 26 23:17:29 2012 +0000
+++ b/mt_resample.m	Mon Feb 27 14:41:24 2012 +0000
@@ -10,7 +10,7 @@
 % mt_init.
 
 function Sys=mt_resample(Sys,K)
-	figure(Sys.fig); 
+	figure(Sys.fig+(K-1)); 
 	cla;
 
 

--- a/private/mc_global_info1.m	Sun Feb 26 23:17:29 2012 +0000
+++ b/private/mc_global_info1.m	Mon Feb 27 14:41:24 2012 +0000
@@ -13,7 +13,7 @@
 function [Info,ergodic]=mc_global_info1(T,tol)
 
 	n=size(T,1);
-	pz=mc_fixpt(T,tol); 
+	pz=mc_fixpt2(T,tol); 
 	ergodic=(size(pz,2)==1);
 	if ergodic && all(isfinite(pz))
 
@@ -60,24 +60,19 @@
 	y=log(x);
 end
 
+function p=mc_fixpt2(Pt,tol)
+	[U,S,V]=svd(Pt-eye(size(Pt)),0);
+	p=V(:,diag(S)<=tol);
+	p=max(0,p*diag(1./sum(p)));
+end
+
+
 % mc_fixpt - fixed point(s) of Markov chain (stationary state distribution)
 %
 % mc_fixpt :: [[K,K]]~'transition matrix'-> [[K,L]].
 function p=mc_fixpt(Pt,tol)
-	% use eigenvalues of transition matrix.
 	[V,D]=eig(Pt);
-
-	% find eigenvalues near 1 to within tolerance
-	k=find(abs(diag(D)-1)<tol);
-	
-	Vk=V(:,k);
-	%sumk=sum(V(:,k));
-	%Vk=V(:,k)*diag(1./sumk); % normalise to unit sum
-	%poss=abs(sumk)>0; % find those with non-zero sum
-	%Vk=V(:,k(poss))*diag(1./sumk(poss)); % normalise to unit sum
-	%p=Vk(:,all(Vk >= -eps)); % eigs with all positive components to within eps.
-
-	%p=max(0,real(Vk*diag(1./sum(Vk))));
-	p=real(Vk*diag(1./sum(Vk)));
-	p(abs(p)<tol/10)=0;
+	p=V(:,abs(diag(D)-1)<tol);
+	p=real(p*diag(1./sum(p)));
+	p(abs(p)<eps)=0;
 end

--- a/private/randnat.m	Sun Feb 26 23:17:29 2012 +0000
+++ b/private/randnat.m	Mon Feb 27 14:41:24 2012 +0000
@@ -1,18 +1,13 @@
-function X=randnat(M,varargin)
-% randnat - Draw Random natural numbers from uniform distribution
+function X=randnat(M)
+% randnat - Draw random natural number from uniform distribution
 %
 % randnat :: 
 %    M:natural   ~'maximum value', 
-%    D:[[1,E]]   ~'size of E-dimensional array'
 % -> [[D]->[M]] ~'size D array of numbers in 1..M'.
 %
 % The size is specified as in the RAND function, except that if
 % only one dimension is given, we create a vector, not a square
 % matrix, eg
-%
-% randnat(10,3) :: [[3,3]->[10]].
-% randnat(20,2,5,7) :: [[2,5,7]->[20]].
-% randnat(20,5) :: [[5,1]->[20]].
 
-X=1+floor(M*rand(tosize(varargin{:})));
+X=1+floor(M*rand);