changeset 0:be936975f254

Initial check in.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mt_calibrate.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,42 @@
+% mt_calibrate - Set spatial calibration of melody triangle system.
+%
+% mt_calibrate ::
+%    mt_system       ~'system to calibrate',
+%    I:[[1,E]->[3]]  ~'list triangle corners being updated',
+%    P:[[3,E]]       ~'X-Y coordinates of the E corners being updated'
+% -> mt_system       ~'calibrated system'.
+
+function Sys=mt_calibrate(Sys,I,PI)
+	
+	if nargin>2, Sys.refpoints(:,I) = PI; end
+	P  = Sys.refpoints;
+	p0 = P(:,1);
+	Sys.map = info_map_fn(p0,inv(P(:,2:3)-repmat(p0,1,2)));
+end
+
+function f=info_map_fn(p0,M)
+	f=@(p)clip_tri(M*(p-p0));
+end
+
+function x=clip_tri(x)
+	if x(2)<0
+		x(2)=0;
+		if x(1)<0, x(1)=0;
+		elseif x(1)>1, x(1)=1;
+		end
+	elseif x(1)<0
+		x(1)=0;
+		if x(2)<0, x(2)=1;
+		elseif x(2)>1, x(2)=1;
+		end
+	else
+		d = [1,1]*x - 1;
+		if d>0
+			x = x - d/2;
+			if x(1)<0, x=[0;1];
+			elseif x(2)<0, x=[1;0];
+			end
+		end
+	end
+end
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mt_ensure.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,17 @@
+% mt_ensure - Ensure that Melody Triange system includes transmats of given size
+%
+% mt_ensure ::
+%    mt_system   ~'melody triangle system',
+%    K:natural   ~'desired size of transition matrices'
+% -> action mt_system.
+%
+% If the system already contains transition matrices of size K, then it is
+% returned unchanged. Otherwise, a set of transition matrices is sampled using
+% the parameters supplied to mt_init. If sampling is done, a scatter plot
+% will be generated as a side effect.
+
+function Sys=mt_ensure(Sys,K)
+	if length(Sys.transmats)<K || isempty(Sys.transmats{K})
+		Sys=mt_resample(Sys,K);
+	end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mt_get_transmat_at.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,37 @@
+% mt_get_transmat_at - Get transition matrix near given point in information space.
+%
+% mt_get_transmat_at ::
+%    mt_system ~'Melody Triangle system structure',
+%    natural   ~'voice Id (1 or greater)',
+%    K:natural ~'size of transition matrix',
+%    real      ~'X-coordinate',
+%    real      ~'Y-coordinate'
+% -> [[K,K]]   ~'transition matrix',
+%    [[1,3]]   ~'information coordinates'.
+%
+% The selected transition matrix will be shown using mt_show_transmat,
+% and it's information coordinates printed in the title of the main scatter plot.
+
+function [T,I]=mt_get_transmat_at(Sys,Id,K,X,Y)
+	normpos = Sys.map([X;Y])';
+	target = normpos *log(K);
+	tmats = Sys.transmats{K};
+	info = Sys.info{K};
+	L = size(info,1);
+
+	% distance from target
+	d2 = sum((info(:,1:2) - repmat(target,L,1)).^2,2);
+	Mask = d2<=max(min(d2)*2,0.01);
+
+	M = find(Mask);
+%	J = M(argmax(info(M,3)));
+	J = M(1+floor(length(M)*rand));
+	T = tmats(:,:,J);
+	I = info(J,:);
+	%Mask=2*int16(Mask); Mask(J)=1;
+	%set(Sys.hScat{K},'CDATA',Mask); % scatc(info,Mask,16); rotate3d on;
+	figure(Sys.fig); title(sprintf('info: %.2f, %.2f, %.2f',I(1),I(2),I(3)));
+	%figure(50); title(sprintf('pos: %.2f, %.2f',normpos(1), normpos(2)));
+	%fprintf('normpos: %.2f, %.2f',normpos(1), normpos(2));
+	mt_show_transmat(Id,T);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mt_init.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,27 @@
+% mt_init - Initialise Melody Triangle system.
+%
+% mt_init :: 
+%    A0:nonneg  ~'parameter for sampling',
+%    B0:nonneg  ~'parameter for sampling',
+%    A1:nonneg  ~'parameter for sampling',
+%    B0:nonneg  ~'parameter for sampling',
+%    L:natural  ~'number of transmats to sample'
+% -> mt_system.
+%
+% Initial system contains no transition matrices - you
+% must call mt_ensure with a particular value of K
+% to sample a set of L transition matrices of that size.
+%
+% Initial calibration is equivalent to:
+%    sys=mt_calibrate(sys, 1:3, [0,1,0;0,0,1]);
+%
+% The figure for scatter plots is fixed to figure 50 for now.
+
+function Sys=mt_init(A0,B0,A1,B1,L)
+	Sys.sample_transmats = @(k)sample_transmat_hdp(A0,B0,A1,B1,k,L);
+	Sys.transmats = {};
+	Sys.info      = {};
+	Sys.refpoints = [0,0;1,0;0,1]';
+	Sys.fig       = 50;
+	Sys = mt_calibrate(Sys);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mt_resample.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,28 @@
+% mt_resample - Sample or resample melody triangle transition matrices of given size
+%
+% mt_resample ::
+%    mt_system  ~'initial system'
+%    natural    ~'size of transition matrices to resample'
+% -> action mt_system.
+%
+% A new set of transition matrices will be sampled and a 3D information space
+% scatter plot generated in the figure determined by the initial call to
+% mt_init.
+
+function Sys=mt_resample(Sys,K)
+	[TT, II] = Sys.sample_transmats(K);
+	figure(Sys.fig); % 'name','info space'); 
+	Sys.transmats{K} = TT;
+	Sys.info{K}      = II;
+	Sys.hScat{K}     = scatc(II,II(:,3),16); 
+	axis on; box on; grid off; 
+	lc=[0.4,0.4,0.4];
+	set(gca,'Color','none');
+	set(gca,'XColor',lc);
+	set(gca,'YColor',lc);
+	set(gca,'ZColor',lc);
+	xlabel('entropy rate');
+	ylabel('redundancy');
+	zlabel('pred-info rate');
+	rotate3d on;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mt_show_transmat.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,11 @@
+% mt_show_transmat - Display image of transition matrix for voice
+%
+% mt_show_transmat ::
+%    natural ~'voice id'
+%    [[K,K]] ~'transition matrix'
+% -> action void.
+
+function mt_show_transmat(Id,T)
+	figure(Id); 
+	imagesc(T,[0,1]); 
+	axis xy;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/mc_fixpt.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,34 @@
+function p=mc_fixpt(Pt)
+% mc_fixpt - fixed point of Markov chain (stationary state distribution)
+%
+% mc_fixpt :: [[K,K]]~'transition matrix'-> [[K]].
+
+
+% for ergodic HMM we want to get stationary distribution 
+% use eigenvalues of transition matrix.
+[V,D]=eig(Pt);
+
+% one of the eigenvalues should be 1
+[dummy,k]=max(-abs(diag(D)-1));
+
+% check that it is really close
+if abs(D(k,k)-1)>0.001
+	disp('mc_fixpt: failed to find eigenvalue of 1 in ');
+	disp(Pt)
+	error('mc_fixpt:noeig','failed to find eigenvalue of 1');
+end
+
+% get eigenvector and flip it over if all negative
+v=V(:,k); if sum(v)<0, v=-v; end
+
+if ~isreal(v)
+	disp('mc_fixpt: discarding complex parts of eigenvector'); 
+	v=real(v);
+end
+if any(v<0)
+	%fprintf('mc_fixpt: discarding negative parts of eigenvector: %s\n',mat2str(v(v<0))); 
+	v=max(0,v);
+end
+
+% final normalisation
+p=v/sum(v);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/mc_global_info.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,55 @@
+function Info=mc_global_info(T,pz)
+% mc_global_info - Three global information measures about stationary MC
+%
+% mc_global_info :: [[N,N]] ~'transmat' -> [[3]].
+%
+% mc_global_info :: 
+%    [[N,N]] ~'transmat',
+%    [[N]] ~'stationary distribution if already known'
+% -> [[3]].
+%
+% The three measures in order are
+%    Entropy rate H(X|Z)
+%    Redundancy I(X,Z)
+%    Predictive information rate I(X,Y|Z)
+
+
+n=size(T,1);
+if nargin<2, pz=mc_fixpt(T); end
+TlogT=T.*slog(T);
+HXZ=-sum(TlogT,1)*pz;
+
+pxz=T; % p(x|z)
+pyxz=repmat(T,[1,1,n]).*repmat(shiftdim(T,-1),[n,1,1]); % p(y,x|z)
+pyz=sum(pyxz,2); % p(y|z)
+pyz(pyz==0)=realmin;
+pxyz=permute(pyxz./repmat(pyz,[1,n,1]),[2,1,3]); % p(x|y,z)
+
+HXYz=-shiftdim(sum(sum(permute(pyxz,[2,1,3]).*slog(pxyz),1),2),2); % H(X|Y,z)
+IXYZ=(-HXYz' - sum(pxz.*slog(pxz)))*pz;
+
+IXZ=max(0,(sum(TlogT)-slog(pz)')*pz); 
+
+
+Info=[HXZ;IXZ;IXYZ];
+if any(~isreal(Info))
+	disp('unreal information');
+	keyboard
+elseif any(Info<0)
+	if any(Info<-eps)
+		disp('negative information');
+		disp(Info);
+		keyboard
+	else
+		Info=max(0,Info);
+	end
+end
+
+function y=slog(x)
+% slog - safe log, x<=0 => slog(x)=0
+%
+% slog :: real -> real.
+
+	x(x<=0)=1;
+	y=log(x);
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/sample_transmat_hdp.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,58 @@
+function [T,I]=sample_transmat_hdp(A0,B0,A1,B1,K,L,T0)
+% sample_transmat_hdp - Sample 1st order Markov transition matrix using HDP
+%
+% sample_transmat_hdp ::
+%    nonneg ~'Shape parameter for base Dirichlet Gamma prior',
+%    nonneg ~'Scale parameter for base Dirichlet Gamma prior',
+%    nonneg ~'Shape parameter for sub Dirichlet Gamma prior',
+%    nonneg ~'Scale parameter for sub Dirichlet Gamma prior',
+%    natural ~'size of alphabet',
+%    natural ~'num of transition matrices to produce',
+% -> [[K,K,L]] ~ 'lots of transition matrices',
+%    [[L,3]]   ~ 'entropy rate, redundancy, and predinfo rate for all'.
+%
+% An extra optional argument of type [[K,K]] can be added, which
+% is used as a sort of bias for the sampled matrices.
+
+	if nargin<7, T0=ones(K,K); end
+
+	Alpha0=gamrnd(A0,B0,1,1,L);
+	Alpha1=gamrnd(A1,B1,1,1,L);
+
+
+	% base distributions
+	H0=sample_dir(repmat(Alpha0/K,[K,1,1]));
+	HT=stoch(repmat(H0,[1,K,1])).*repmat(T0/sum(T0(:)),[1,1,L]);
+	T=sample_dir(repmat(Alpha1,[K,K,1]).*HT);
+
+	if any(sum(H0)<0.99)
+		disp('normalisation problem in H')
+		keyboard
+	end
+	if nargout>=2
+		I=zeros(K,3);
+		for i=1:L
+			I(i,:)=mc_global_info(T(:,:,i));
+		end
+	end
+end
+
+% taken from stats/models/dirichlet.m>stoch_gamma
+function y=sample_dir(beta)
+	x=randg(beta);
+	z=sum(x,1);
+	y=x./repmat(z,size(x,1),1);
+	bad=find(z==0);
+	for i=1:length(bad)
+		beta_i=beta(:,bad(i)); % misbehaving params
+		k=(beta_i==max(beta_i)); % entries with max beta
+		y(:,bad(i))=mnrnd(1,stoch(k));
+	end
+end
+
+function [Y,Z]=stoch(X),
+	Z=sum(X,1);
+	Y=X./repmat(Z,size(X,1),1);
+end
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/scatc.m	Wed Feb 01 14:06:37 2012 +0000
@@ -0,0 +1,43 @@
+function h=scatc(x,C,S,varargin)
+% scatc - 2 or 3D scatter plot with colours and sizes
+%
+% scatc :: [[N,E]] ~'N points in E space', [[N]]~'colours' -> handle.
+% scatc :: 
+%    [[N,E]] ~'N points in E space', 
+%    [[N]]   ~'colours', 
+%    [[N]] | real  ~'array of sizes or single marker size'
+% -> handle.
+%
+% If E<3, does 2D scatter, otherwise, does 3D scatter plot.
+% Draws filled circles. 
+
+if size(x,2)<=4 && size(x,1)>4, x=x'; end
+if nargin<3 || isempty(S),
+	% use default marker size
+	S=get(gca,'DefaultLineMarkerSize').^2; 
+elseif length(S)==1,
+	S=S*ones(size(x,2),1);
+end
+
+if size(x,1)<3,
+	h=scatter(x(1,:),x(2,:),S,C);
+else
+	h=scatter3(x(1,:),x(2,:),x(3,:),S,C);
+end
+
+set(h, ...
+	'MarkerFaceColor','flat', ...
+	'MarkerEdgeColor','k', ...
+	'LineWidth',0.2, ...
+	varargin{:});
+
+axis equal;
+if size(x,1)>2
+	axis vis3d;
+	axis off;
+	set(gca,'DrawMode','fast');
+	set(gca,'Projection','perspective');
+else
+	box on;
+	% axis off;
+end