wolffd@0: function bnet = mk_hhmm3(varargin)
wolffd@0: % MK_HHMM3 Make a 3 level Hierarchical HMM
wolffd@0: % bnet = mk_hhmm3(...)
wolffd@0: %
wolffd@0: % 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs.
wolffd@0: % This enforces sub-models (which differ only in their Q1 index) to be shared.
wolffd@0: % Also, we enforce the fact that each model always starts in its initial state
wolffd@0: % and only finishes in its final state. However, the prob. of finishing (as opposed to
wolffd@0: % self-transitioning to the final state) can be learned.
wolffd@0: % The fact that we always finish from the same state means we do not need to condition
wolffd@0: % F(i) on Q(i-1), since finishing prob is indep of calling context.
wolffd@0: %
wolffd@0: % The DBN is the same as Fig 10 in my tech report.
wolffd@0: %
wolffd@0: %   Q1 ---------->  Q1
wolffd@0: %   |              / |
wolffd@0: %   |             /  |
wolffd@0: %   |  F2 -------    |
wolffd@0: %   |  ^         \   |
wolffd@0: %   | /|          \  |
wolffd@0: %   v  |           v v
wolffd@0: %   Q2-| -------->   Q2
wolffd@0: %  /|  |             ^
wolffd@0: % / |  |            /|
wolffd@0: % | |  F3 ---------/ |
wolffd@0: % | |  ^           \ |
wolffd@0: % | v /              v
wolffd@0: % | Q3 ----------->  Q3
wolffd@0: % |  |    
wolffd@0: % \  | 
wolffd@0: %  v v    
wolffd@0: %   O
wolffd@0: %
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: % Q3 (slice 1) is clamped to always start in state 1.
wolffd@0: % F3 by default will only finish if Q3 is in its last state (F3 is a tabular_CPD)
wolffd@0: % F2 by default gets the default hhmmF_CPD params.
wolffd@0: % Q1:Q3 (slice 2) by default gets the default hhmmQ_CPD params.
wolffd@0: % O by default gets the default tabular/Gaussian params.
wolffd@0: %
wolffd@0: % Optional arguments in name/value format [default]
wolffd@0: %
wolffd@0: % Qsizes      - sizes at each level [ none ]
wolffd@0: % Osize       - size of O node [ none ]
wolffd@0: % discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
wolffd@0: % Oargs       - cell array of args to pass to the O CPD  [ {} ]
wolffd@0: % Q1args      - args to be passed to constructor for Q1 (slice 2) [ {} ]
wolffd@0: % Q2args      - args to be passed to constructor for Q2 (slice 2) [ {} ]
wolffd@0: % Q3args      - args to be passed to constructor for Q3 (slice 2) [ {} ]
wolffd@0: % F2args       - args to be passed to constructor for F2 [ {} ]
wolffd@0: % F3args       - args to be passed to constructor for F3 [ {'CPT', finish in last Q3 state} ]
wolffd@0: %
wolffd@0: 
wolffd@0: ss = 6; D = 3;
wolffd@0: Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
wolffd@0: Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
wolffd@0: names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'};
wolffd@0: 
wolffd@0: intra = zeros(ss);
wolffd@0: intra(Q1, Q2) = 1;
wolffd@0: intra(Q2, [F2 Q3 obs]) = 1;
wolffd@0: intra(Q3, [F3 obs]) = 1;
wolffd@0: intra(F3, F2) = 1;
wolffd@0: 
wolffd@0: inter = zeros(ss);
wolffd@0: inter(Q1,Q1) = 1;
wolffd@0: inter(Q2,Q2) = 1;
wolffd@0: inter(Q3,Q3) = 1;
wolffd@0: inter(F2,[Q1 Q2]) = 1;
wolffd@0: inter(F3,[Q2 Q3]) = 1;
wolffd@0: 
wolffd@0: 
wolffd@0: % get sizes of nodes
wolffd@0: args = varargin;
wolffd@0: nargs = length(args);
wolffd@0: Qsizes = [];
wolffd@0: Osize = 0;
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'Qsizes', Qsizes = args{i+1}; 
wolffd@0:    case 'Osize', Osize = args{i+1}; 
wolffd@0:   end
wolffd@0: end
wolffd@0: if isempty(Qsizes), error('must specify Qsizes'); end
wolffd@0: if Osize==0, error('must specify Osize'); end
wolffd@0:   
wolffd@0: % set default params
wolffd@0: discrete_obs = 0;
wolffd@0: Oargs = {};
wolffd@0: Q1args = {};
wolffd@0: Q2args = {};
wolffd@0: Q3args = {};
wolffd@0: F2args = {};
wolffd@0: 
wolffd@0: % P(Q3, F3)
wolffd@0: CPT = zeros(Qsizes(3), 2);
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(Qsizes(3), 2) = p; % except the last one
wolffd@0: CPT(Qsizes(3), 1) = 1-p;
wolffd@0: F3args = {'CPT', CPT};
wolffd@0: 
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'discrete_obs', discrete_obs = args{i+1}; 
wolffd@0:    case 'Oargs',        Oargs = args{i+1};
wolffd@0:    case 'Q1args',       Q1args = args{i+1};
wolffd@0:    case 'Q2args',       Q2args = args{i+1};
wolffd@0:    case 'Q3args',       Q3args = args{i+1};
wolffd@0:    case 'F2args',       F2args = args{i+1};
wolffd@0:    case 'F3args',       F3args = args{i+1};
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: ns = zeros(1,ss);
wolffd@0: ns(Qnodes) = Qsizes;
wolffd@0: ns(obs) = Osize;
wolffd@0: ns(Fnodes) = 2;
wolffd@0: 
wolffd@0: dnodes = [Qnodes Fnodes];
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,ns(Q1)));
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: % each model always starts in state 1
wolffd@0: CPT = zeros(ns(Q2), ns(Q3));
wolffd@0: CPT(:, 1) = 1.0;
wolffd@0: bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', CPT, 'adjustable', 0);
wolffd@0: 
wolffd@0: bnet.CPD{eclass(F2,1)}  = hhmmF_CPD(bnet, F2, Qnodes, 2, D, F2args{:});
wolffd@0: 
wolffd@0: bnet.CPD{eclass(F3,1)}  = tabular_CPD(bnet, F3, F3args{:});
wolffd@0: 
wolffd@0: if discrete_obs
wolffd@0:   bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:});
wolffd@0: else
wolffd@0:   bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:});
wolffd@0: end
wolffd@0: 
wolffd@0: % SLICE 2
wolffd@0: 
wolffd@0: bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, Q1args{:});
wolffd@0: bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:});
wolffd@0: bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, Q3args{:});