wolffd@0: 
wolffd@0: seed = 0;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: 
wolffd@0: discrete_obs = 1;
wolffd@0: topright = 0;
wolffd@0: 
wolffd@0: Qsizes = [2 4 2];
wolffd@0: D = 3;
wolffd@0: Qnodes = 1:D;
wolffd@0: startprob = cell(1,D);
wolffd@0: transprob = cell(1,D);
wolffd@0: termprob = cell(1,D);
wolffd@0: 
wolffd@0: % LEVEL 1
wolffd@0: 
wolffd@0: startprob{1} = 'ergodic';
wolffd@0: transprob{1} = 'ergodic';
wolffd@0: 
wolffd@0: % LEVEL 2
wolffd@0: 
wolffd@0: startprob{2} = zeros(2, 4);
wolffd@0: startprob{2}(1, :) = [1 0 0 0];
wolffd@0: if topright
wolffd@0:   startprob{2}(2, :) = [0 0 1 0];
wolffd@0: else
wolffd@0:   startprob{2}(2, :) = [0 1 0 0];
wolffd@0: end
wolffd@0: 
wolffd@0: transprob{2} = zeros(4, 2, 4);
wolffd@0: 
wolffd@0: transprob{2}(:,1,:) = [0 1 0 0
wolffd@0: 		       0 0 1 0
wolffd@0: 		       0 0 0 1
wolffd@0: 		       0 0 0 1]; % 4->e
wolffd@0: if topright
wolffd@0:   transprob{2}(:,2,:) = [0 0 0 1
wolffd@0: 		    1 0 0 0
wolffd@0: 		    0 1 0 0
wolffd@0: 		    0 0 0 1]; % 4->e
wolffd@0: else
wolffd@0:   transprob{2}(:,2,:) = [0 0 0 1
wolffd@0: 		    1 0 0 0
wolffd@0: 		    0 0 1 0 % 3->e
wolffd@0: 		    0 0 1 0];
wolffd@0: end
wolffd@0: 
wolffd@0: %termprob{2} = 'rightstop';
wolffd@0: termprob{2} = zeros(2,4,2);
wolffd@0: pfin = 0.8;
wolffd@0: termprob{2}(1,:,2) = [0 0 0 pfin]; % finish in state 4 (DU)
wolffd@0: termprob{2}(1,:,1) = 1 - [0 0 0 pfin];
wolffd@0: if topright
wolffd@0:   termprob{2}(2,:,2) = [0 0 0 pfin];
wolffd@0:   termprob{2}(2,:,1) = 1 - [0 0 0 pfin];
wolffd@0: else
wolffd@0:   termprob{2}(2,:,2) = [0 0 pfin 0];  % finish in state 3 (RL)
wolffd@0:   termprob{2}(2,:,1) = 1 - [0 0 pfin 0];
wolffd@0: end
wolffd@0: 
wolffd@0: % LEVEL 3
wolffd@0: 
wolffd@0: startprob{3} = 'leftstart';
wolffd@0: transprob{3}  = 'leftright';
wolffd@0: termprob{3} = 'rightstop';
wolffd@0: 
wolffd@0: 
wolffd@0: % OBS LEVEl
wolffd@0: 
wolffd@0: if discrete_obs
wolffd@0:   chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
wolffd@0:   L=find(chars=='L'); l=find(chars=='l');
wolffd@0:   U=find(chars=='U'); u=find(chars=='u');
wolffd@0:   R=find(chars=='R'); r=find(chars=='r');
wolffd@0:   D=find(chars=='D'); d=find(chars=='d');
wolffd@0:   Osize = length(chars);
wolffd@0:   
wolffd@0:   obsprob = zeros([4 2 Osize]);
wolffd@0:   %       Q2 Q3 O
wolffd@0:   obsprob(1, 1, L) =  1.0;
wolffd@0:   obsprob(1, 2, l) =  1.0;
wolffd@0:   obsprob(2, 1, U) =  1.0;
wolffd@0:   obsprob(2, 2, u) =  1.0;
wolffd@0:   obsprob(3, 1, R) =  1.0;
wolffd@0:   obsprob(3, 2, r) =  1.0;
wolffd@0:   obsprob(4, 1, D) =  1.0;
wolffd@0:   obsprob(4, 2, d) =  1.0;
wolffd@0:   
wolffd@0:   Oargs = {'CPT', obsprob};
wolffd@0: else
wolffd@0:   Osize = 2;
wolffd@0:   mu = zeros(2, 4, 2);
wolffd@0:   noise = 0;
wolffd@0:   scale = 10;
wolffd@0:   for q3=1:2
wolffd@0:     mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1);
wolffd@0:   end
wolffd@0:   for q3=1:2
wolffd@0:     mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1);
wolffd@0:   end
wolffd@0:   for q3=1:2
wolffd@0:     mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1);
wolffd@0:   end
wolffd@0:   for q3=1:2
wolffd@0:     mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1);
wolffd@0:   end
wolffd@0:   Sigma = repmat(reshape(0.01*eye(2), [2 2 1 1 ]), [1 1 4 2]);
wolffd@0:   Oargs = {'mean', mu, 'cov', Sigma};
wolffd@0: end
wolffd@0: 
wolffd@0: bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ...
wolffd@0: 	       'Oargs', Oargs, 'Ops', Qnodes(2:3), ...
wolffd@0: 	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
wolffd@0: 
wolffd@0: if discrete_obs
wolffd@0:   Tmax = 30;
wolffd@0: else
wolffd@0:   Tmax = 200;
wolffd@0: end
wolffd@0: usecell = ~discrete_obs;
wolffd@0: Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
wolffd@0: Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
wolffd@0: 
wolffd@0: for seqi=1:3
wolffd@0:   evidence = sample_dbn(bnet, Tmax, usecell, 'stop_sampling_F2');      
wolffd@0:   T = size(evidence, 2)
wolffd@0:   if discrete_obs
wolffd@0:     pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars);
wolffd@0:   else
wolffd@0:     pos = zeros(2,T+1);
wolffd@0:     delta = cell2num(evidence(Onode,:));
wolffd@0:     clf
wolffd@0:     hold on
wolffd@0:     cols = {'r', 'g', 'k', 'b'};
wolffd@0:     boundary = cell2num(evidence(F3,:))-1;
wolffd@0:     coli = 1;
wolffd@0:     for t=2:T+1
wolffd@0:       pos(:,t) = pos(:,t-1) + delta(:,t-1);
wolffd@0:       plot(pos(1,t), pos(2,t), sprintf('%c.', cols{coli}));
wolffd@0:       if boundary(t-1)
wolffd@0: 	coli = coli + 1;
wolffd@0: 	coli = mod(coli-1, length(cols)) + 1;
wolffd@0:       end
wolffd@0:     end
wolffd@0:     %plot(pos(1,:), pos(2,:), '.')
wolffd@0:     %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, []);
wolffd@0:     pause
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: eclass = bnet.equiv_class;
wolffd@0: S=struct(bnet.CPD{eclass(Q2,2)});
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: