wolffd@0: % Try to learn a 3 level HHMM similar to mk_square_hhmm
wolffd@0: % from hand-drawn squares.
wolffd@0: 
wolffd@0: % Because startprob should be shared for t=1:T,
wolffd@0: % but in the DBN is shared for t=2:T, we train using a single long sequence.
wolffd@0: 
wolffd@0: discrete_obs = 0;
wolffd@0: supervised = 1;
wolffd@0: obs_finalF2 = 0;
wolffd@0: % It is not possible to observe F2 if we learn
wolffd@0: % because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume
wolffd@0: % the F nodes are always hidden (for speed).
wolffd@0: % However, for generating, we might want to set the final F2=true
wolffd@0: % to force all subroutines to finish.
wolffd@0: 
wolffd@0: seed = 1;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: 
wolffd@0: bnet = mk_square_hhmm(discrete_obs, 0);
wolffd@0:  
wolffd@0: ss = 6;
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: Qsizes = [2 4 1];
wolffd@0:   
wolffd@0: if supervised
wolffd@0:   bnet.observed = [Q1 Q2 Onode];
wolffd@0: else
wolffd@0:   bnet.observed = [Onode];
wolffd@0: end
wolffd@0: 
wolffd@0: if obs_finalF2
wolffd@0:   engine = jtree_dbn_inf_engine(bnet);
wolffd@0:   % can't use ndx version because sometimes F2 is hidden, sometimes observed
wolffd@0:   error('can''t observe F when learning')
wolffd@0: else
wolffd@0:   if supervised
wolffd@0:     engine = jtree_ndx_dbn_inf_engine(bnet);
wolffd@0:   else
wolffd@0:     engine = jtree_hmm_inf_engine(bnet);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: load 'square4_cases' % cases{seq}{i,t} for i=1:ss 
wolffd@0: %plot_square_hhmm(cases{1})
wolffd@0: %long_seq = cat(2, cases{:});
wolffd@0: train_cases = cases(1:2);
wolffd@0: long_seq = cat(2, train_cases{:});
wolffd@0: if ~supervised
wolffd@0:   T = size(long_seq,2);
wolffd@0:   for t=1:T
wolffd@0:     long_seq{Q1,t} = [];
wolffd@0:     long_seq{Q2,t} = [];
wolffd@0:   end
wolffd@0: end
wolffd@0: [bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 2);
wolffd@0: 
wolffd@0: eclass = bnet2.equiv_class;
wolffd@0: CPDO=struct(bnet2.CPD{eclass(Onode,1)});
wolffd@0: mu = CPDO.mean;
wolffd@0: Sigma = CPDO.cov;
wolffd@0: CPDO_full = CPDO;
wolffd@0: 
wolffd@0: % force diagonal covs after training
wolffd@0: for k=1:size(Sigma,3)
wolffd@0:   Sigma(:,:,k) = diag(diag(Sigma(:,:,k)));
wolffd@0: end
wolffd@0: bnet2.CPD{6} = set_fields(bnet.CPD{6}, 'cov', Sigma);
wolffd@0: 
wolffd@0: if 0
wolffd@0:   % visualize each model by concatenating means for each model for nsteps in a row
wolffd@0:   nsteps = 5;
wolffd@0:   ev = cell(ss, nsteps*prod(Qsizes(2:3)));
wolffd@0:   t = 1;
wolffd@0:   for q2=1:Qsizes(2)
wolffd@0:     for q3=1:Qsizes(3)
wolffd@0:       for i=1:nsteps
wolffd@0: 	ev{Onode,t} = mu(:,q2,q3);
wolffd@0: 	ev{Q2,t} = q2;
wolffd@0: 	t = t + 1;
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0:   plot_square_hhmm(ev)      
wolffd@0: end
wolffd@0: 
wolffd@0: % bnet3 is the same as the learned model, except we will use it in testing mode
wolffd@0: if supervised
wolffd@0:   bnet3 = bnet2;
wolffd@0:   bnet3.observed = [Onode];
wolffd@0:   engine3 = hmm_inf_engine(bnet3);
wolffd@0:   %engine3 = jtree_ndx_dbn_inf_engine(bnet3);
wolffd@0: else
wolffd@0:   bnet3 = bnet2;
wolffd@0:   engine3 = engine2;
wolffd@0: end
wolffd@0: 
wolffd@0: if 0
wolffd@0:   % segment whole sequence
wolffd@0:   mpe = calc_mpe_dbn(engine3, long_seq);
wolffd@0:   pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []);
wolffd@0: end
wolffd@0: 
wolffd@0: % segment each sequence
wolffd@0: test_cases = cases(3:4);
wolffd@0: for i=1:2
wolffd@0:   ev = test_cases{i};
wolffd@0:   T = size(ev, 2);
wolffd@0:   for t=1:T
wolffd@0:     ev{Q1,t} = [];
wolffd@0:     ev{Q2,t} = [];
wolffd@0:   end
wolffd@0:   %mpe = calc_mpe_dbn(engine3, ev);
wolffd@0:   mpe = find_mpe(engine3, ev)
wolffd@0:   subplot(1,2,i)
wolffd@0:   plot_square_hhmm(mpe)      
wolffd@0:   %pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []);
wolffd@0:   q1s = cell2num(mpe(Q1,:));
wolffd@0:   h = hist(q1s, 1:Qsizes(1));
wolffd@0:   map_q1 = argmax(h);
wolffd@0:   str = sprintf('test seq %d is of type %d\n', i, map_q1);
wolffd@0:   title(str)
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: if 0
wolffd@0: % Estimate gotten by couting transitions in the labelled data
wolffd@0: % Note that a self transition shouldnt count if F2=off.
wolffd@0: Q2ev = cell2num(ev(Q2,:));
wolffd@0: Q2a = Q2ev(1:end-1);
wolffd@0: Q2b = Q2ev(2:end);
wolffd@0: counts = compute_counts([Q2a; Q2b], [4 4]);
wolffd@0: end
wolffd@0: 
wolffd@0: eclass = bnet2.equiv_class;
wolffd@0: CPDQ1=struct(bnet2.CPD{eclass(Q1,2)});
wolffd@0: CPDQ2=struct(bnet2.CPD{eclass(Q2,2)});
wolffd@0: CPDQ3=struct(bnet2.CPD{eclass(Q3,2)});
wolffd@0: CPDF2=struct(bnet2.CPD{eclass(F2,1)});
wolffd@0: CPDF3=struct(bnet2.CPD{eclass(F3,1)});
wolffd@0: 
wolffd@0: 
wolffd@0: A=add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2));
wolffd@0: squeeze(A(:,1,:));
wolffd@0: CPDQ2.startprob;
wolffd@0:  
wolffd@0: if 0
wolffd@0: S=struct(CPDF2.sub_CPD_term);
wolffd@0: S.nsamples
wolffd@0: reshape(S.counts, [2 4 2])
wolffd@0: end