wolffd@0: function bnet = mk_hhmm2(varargin)
wolffd@0: % MK_HHMM2 Make a 2 level Hierarchical HMM
wolffd@0: % bnet = mk_hhmm2(...)
wolffd@0: %
wolffd@0: % 2-layer hierarchical HMM  (node numbers in parens)
wolffd@0: %
wolffd@0: %   Q1(1) ---------> Q1(5)
wolffd@0: % /  | \            / |
wolffd@0: % |  |  v          /  |
wolffd@0: % |  |  F2(3) --- /   |
wolffd@0: % |  |  ^         \   |
wolffd@0: % |  | /           \  |
wolffd@0: % |  v              \ v
wolffd@0: % |  Q2(2)--------> Q2 (6)
wolffd@0: % |  |    
wolffd@0: % \  | 
wolffd@0: %  v v    
wolffd@0: %   O(4)
wolffd@0: %
wolffd@0: %
wolffd@0: % Optional arguments [default]
wolffd@0: %
wolffd@0: % discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
wolffd@0: % obsCPT       - CPT(o,q1,q2) params for O ['rnd']
wolffd@0: % mu           - mu(:,q1,q2) params for O [ [] ]
wolffd@0: % Sigma        - Sigma(:,q1,q2) params for O [ [] ]
wolffd@0: %
wolffd@0: % F2toQ1       - 1 if Q2 is an hhmm_CPD, 0 if F2 -> Q2 arc is absent, so level 2 never resets [1]
wolffd@0: % Q1args        - arguments to be passed to the constructors for Q1(t=2) [ {} ]
wolffd@0: % Q2args        - arguments to be passed to the constructors for Q2(t=2) [ {} ]
wolffd@0: %
wolffd@0: % F2 only turns on (wp 0.5) when Q2 enters its final state.
wolffd@0: % Q1 (slice 1) is clamped to be uniform.
wolffd@0: % Q2 (slice 1) is clamped to always start in state 1.
wolffd@0: 
wolffd@0: [os nmodels nstates] = size(mu);
wolffd@0: 
wolffd@0: ss = 4;
wolffd@0: Q1 = 1; Q2 = 2; F2 = 3; obs = 4;
wolffd@0: Qnodes = [Q1 Q2];
wolffd@0: names = {'Q1', 'Q2', 'F2', 'obs'};
wolffd@0: intra = zeros(ss);
wolffd@0: intra(Q1, [Q2 F2 obs]) = 1;
wolffd@0: intra(Q2, [F2 obs]) = 1;
wolffd@0: 
wolffd@0: inter = zeros(ss);
wolffd@0: inter(Q1,Q1) = 1;
wolffd@0: inter(F2,Q1) = 1;
wolffd@0: if F2toQ2
wolffd@0:   inter(F2,Q2)=1;
wolffd@0: end
wolffd@0: inter(Q2,Q2) = 1;
wolffd@0: 
wolffd@0: ns = zeros(1,ss);
wolffd@0: 
wolffd@0: ns(Q1) = nmodels;
wolffd@0: ns(Q2) = nstates;
wolffd@0: ns(F2) = 2;
wolffd@0: ns(obs) = os;
wolffd@0: 
wolffd@0: dnodes = [Q1 Q2 F2];
wolffd@0: if discrete_obs
wolffd@0:   dnodes = [dnodes obs];
wolffd@0: end
wolffd@0: onodes = [obs];
wolffd@0: 
wolffd@0: bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
wolffd@0: eclass = bnet.equiv_class;
wolffd@0: 
wolffd@0: % SLICE 1
wolffd@0: 
wolffd@0: % We clamp untied nodes in the first slice, since their params can't be estimated
wolffd@0: % from just one sequence
wolffd@0: 
wolffd@0: % uniform prior on initial model
wolffd@0: CPT = normalise(ones(1,nmodels));
wolffd@0: bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', CPT, 'adjustable', 0);
wolffd@0: 
wolffd@0: % each model always starts in state 1
wolffd@0: CPT = zeros(ns(Q1), ns(Q2));
wolffd@0: CPT(:, 1) = 1.0;
wolffd@0: bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', CPT, 'adjustable', 0);
wolffd@0: 
wolffd@0: % Termination probability
wolffd@0: CPT = zeros(ns(Q1), ns(Q2), 2);
wolffd@0: if 1
wolffd@0:   % Each model can only terminate in its final state.
wolffd@0:   % 0 params will remain 0 during EM, thus enforcing this constraint.
wolffd@0:   CPT(:, :, 1) = 1.0; % all states turn F off ...
wolffd@0:   p = 0.5;
wolffd@0:   CPT(:, ns(Q2), 2) = p; % except the last one
wolffd@0:   CPT(:, ns(Q2), 1) = 1-p;
wolffd@0: end
wolffd@0: bnet.CPD{eclass(F2,1)}  = tabular_CPD(bnet, F2, 'CPT', CPT);
wolffd@0: 
wolffd@0: if discrete_obs
wolffd@0:   bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, obs_args{:});
wolffd@0: else
wolffd@0:   bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, obs_args{:});
wolffd@0: end
wolffd@0: 
wolffd@0: % SLICE 2
wolffd@0: 
wolffd@0: 
wolffd@0: bnet.CPD{eclass(Q1,2)} = hhmm_CPD(bnet, Q1+ss, Qnodes, 1, D, 'args', Q1args);
wolffd@0: 
wolffd@0: if F2toQ2
wolffd@0:   bnet.CPD{eclass(Q2,2)} = hhmmQD_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:});
wolffd@0: else
wolffd@0:   bnet.CPD{eclass(Q2,2)} = tabular_CPD(bnet, Q2+ss, Q2args{:});
wolffd@0: end