diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/examples/dynamic/HHMM')
76 files changed, 5878 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries new file mode 100644 index 00000000..6c5b3923 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries @@ -0,0 +1,9 @@ +/abcd_hhmm.m/1.1.1.1/Sat Sep 21 21:37:54 2002// +/add_hhmm_end_state.m/1.1.1.1/Wed May 29 15:59:54 2002// +/hhmm_jtree_clqs.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_hhmm.m/1.1.1.1/Sat Sep 21 20:58:06 2002// +/mk_hhmm_topo.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_hhmm_topo_F1.m/1.1.1.1/Wed May 29 15:59:54 2002// +/pretty_print_hhmm_parse.m/1.1.1.1/Wed May 29 15:59:54 2002// +/remove_hhmm_end_state.m/1.1.1.1/Mon Dec 16 19:16:50 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries.Log new file mode 100644 index 00000000..1b0fe64f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Entries.Log @@ -0,0 +1,5 @@ +A D/Map//// +A D/Mgram//// +A D/Motif//// +A D/Old//// +A D/Square//// diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Repository new file mode 100644 index 00000000..2ba8b9e6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Entries new file mode 100644 index 00000000..dc32f52b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Entries @@ -0,0 +1,6 @@ +/disp_map_hhmm.m/1.1.1.1/Tue Sep 24 22:45:56 2002// +/learn_map.m/1.1.1.1/Sat Jan 11 18:48:46 2003// +/mk_map_hhmm.m/1.1.1.1/Tue Sep 24 10:49:52 2002// +/mk_rnd_map_hhmm.m/1.1.1.1/Tue Sep 24 22:13:48 2002// +/sample_from_map.m/1.1.1.1/Tue Sep 24 13:02:30 2002// +D/Old//// diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Repository new file mode 100644 index 00000000..66b47bbc --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Map diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Entries new file mode 100644 index 00000000..6079d451 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Entries @@ -0,0 +1,2 @@ +/mk_map_hhmm.m/1.1.1.1/Tue Sep 24 07:02:44 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Repository new file mode 100644 index 00000000..354057a9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Map/Old diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/mk_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/mk_map_hhmm.m new file mode 100644 index 00000000..7b646745 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/Old/mk_map_hhmm.m @@ -0,0 +1,156 @@ +function bnet = mk_map_hhmm(varargin) + +% p is the prob of a successful move (defines the reliability of motors) +p = 1; +num_obs_nodes = 1; + +for i=1:2:length(varargin) + switch varargin{i}, + case 'p', p = varargin{i+1}; + case 'numobs', num_obs_node = varargin{i+1}; + end +end + + +q = 1-p; + +% assign numbers to the nodes in topological order +U = 1; A = 2; C = 3; F = 4; O = 5; + +% create graph structure + +ss = 5; % slice size +intra = zeros(ss,ss); +intra(U,F)=1; +intra(A,[C F O])=1; +intra(C,[F O])=1; + +inter = zeros(ss,ss); +inter(U,[A C])=1; +inter(A,[A C])=1; +inter(F,[A C])=1; +inter(C,C)=1; + +% node sizes +ns = zeros(1,ss); +ns(U) = 2; % left/right +ns(A) = 2; +ns(C) = 3; +ns(F) = 2; +ns(O) = 5; % we will assign each state a unique symbol +l = 1; r = 2; % left/right +L = 1; R = 2; + +% Make the DBN +bnet = mk_dbn(intra, inter, ns, 'observed', O); +eclass = bnet.equiv_class; + + + +% Define CPDs for slice 1 +% We clamp all of them, i.e., do not try to learn them. + +% uniform probs over actions (the input could be chosen from a policy) +bnet.CPD{eclass(U,1)} = tabular_CPD(bnet, U, 'CPT', mk_stochastic(ones(ns(U),1)), ... + 'adjustable', 0); + +% uniform probs over starting abstract state +bnet.CPD{eclass(A,1)} = tabular_CPD(bnet, A, 'CPT', mk_stochastic(ones(ns(A),1)), ... + 'adjustable', 0); + +% Uniform probs over starting concrete state, modulo the fact +% that corridor 2 is only of length 2. +CPT = zeros(ns(A), ns(C)); % CPT(i,j) = P(C starts in j | A=i) +CPT(1, :) = [1/3 1/3 1/3]; +CPT(2, :) = [1/2 1/2 0]; +bnet.CPD{eclass(C,1)} = tabular_CPD(bnet, C, 'CPT', CPT, 'adjustable', 0); + +% Termination probs +CPT = zeros(ns(U), ns(A), ns(C), ns(F)); +CPT(r,1,1,:) = [1 0]; +CPT(r,1,2,:) = [1 0]; +CPT(r,1,3,:) = [q p]; +CPT(r,2,1,:) = [1 0]; +CPT(r,2,2,:) = [q p]; +CPT(l,1,1,:) = [q p]; +CPT(l,1,2,:) = [1 0]; +CPT(l,1,3,:) = [1 0]; +CPT(l,2,1,:) = [q p]; +CPT(l,2,2,:) = [1 0]; + +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT); + + +% Assign each state a unique observation +CPT = zeros(ns(A), ns(C), ns(O)); +CPT(1,1,1)=1; +CPT(1,2,2)=1; +CPT(1,3,3)=1; +CPT(2,1,4)=1; +CPT(2,2,5)=1; +%CPT(2,3,:) undefined + +bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT); + + +% Define the CPDs for slice 2 + +% Abstract + +% Since the top level never resets, the starting distribution is irrelevant: +% A2 will be determined by sampling from transmat(A1,:). +% But the code requires we specify it anyway; we make it all 0s, a dummy value. +startprob = zeros(ns(U), ns(A)); + +transmat = zeros(ns(U), ns(A), ns(A)); +transmat(R,1,:) = [q p]; +transmat(R,2,:) = [0 1]; +transmat(L,1,:) = [1 0]; +transmat(L,2,:) = [p q]; + +% Qps are the parents we condition the parameters on, in this case just +% the past action. +bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ... + 'startprob', startprob, 'transprob', transmat); + + + +% Concrete + +transmat = zeros(ns(C), ns(U), ns(A), ns(C)); +transmat(1,r,1,:) = [q p 0.0]; +transmat(2,r,1,:) = [0.0 q p]; +transmat(3,r,1,:) = [0.0 0.0 1.0]; +transmat(1,r,2,:) = [q p 0.0]; +transmat(2,r,2,:) = [0.0 1.0 0.0]; +% +transmat(1,l,1,:) = [1.0 0.0 0.0]; +transmat(2,l,1,:) = [p q 0.0]; +transmat(3,l,1,:) = [0.0 p q]; +transmat(1,l,2,:) = [1.0 0.0 0.0]; +transmat(2,l,2,:) = [p q 0.0]; + +% Add a new dimension for A(t-1), by copying old vals, +% so the matrix is the same size as startprob + + +transmat = reshape(transmat, [ns(C) ns(U) ns(A) 1 ns(C)]); +transmat = repmat(transmat, [1 1 1 ns(A) 1]); + +% startprob(C(t-1), U(t-1), A(t-1), A(t), C(t)) +startprob = zeros(ns(C), ns(U), ns(A), ns(A), ns(C)); +startprob(1,L,1,1,:) = [1.0 0.0 0.0]; +startprob(3,R,1,2,:) = [1.0 0.0 0.0]; +startprob(3,R,1,1,:) = [0.0 0.0 1.0]; +% +startprob(1,L,2,1,:) = [0.0 0.0 010]; +startprob(2,L,2,1,:) = [1.0 0.0 0.0]; +startprob(2,R,2,2,:) = [0.0 1.0 0.0]; + +% want transmat(U,A,C,At,Ct), ie. in topo order +transmat = permute(transmat, [2 3 1 4 5]); +startprob = permute(startprob, [2 3 1 4 5]); +bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ... + 'startprob', startprob, 'transprob', transmat); + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/disp_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/disp_map_hhmm.m new file mode 100644 index 00000000..0aadf2bb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/disp_map_hhmm.m @@ -0,0 +1,13 @@ +function disp_map_hhmm(bnet) + +eclass = bnet.equiv_class; +U = 1; A = 2; C = 3; F = 4; + +S = struct(bnet.CPD{eclass(A,2)}); +disp('abstract trans') +dispcpt(S.transprob) + +S = struct(bnet.CPD{eclass(C,2)}); +disp('concrete trans for go left') % UAC AC +dispcpt(squeeze(S.transprob(1,:,:,:,:))) + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/learn_map.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/learn_map.m new file mode 100644 index 00000000..ac36586a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/learn_map.m @@ -0,0 +1,40 @@ +seed = 1; +rand('state', seed); +randn('state', seed); + +obs_model = 'unique'; % each cell has a unique label (essentially fully observable) +%obs_model = 'four'; % each cell generates 4 observations, NESW + +% Generate the true network, and a randomization of it +realnet = mk_map_hhmm('p', 0.9, 'obs_model', obs_model); +rndnet = mk_rnd_map_hhmm('obs_model', obs_model); +eclass = realnet.equiv_class; +U = 1; A = 2; C = 3; F = 4; onodes = 5; + +ss = realnet.nnodes_per_slice; +T = 100; +evidence = sample_dbn(realnet, 'length', T); +ev = cell(ss,T); +ev(onodes,:) = evidence(onodes,:); + +infeng = jtree_dbn_inf_engine(rndnet); + +if 0 +% suppose we do not observe the final finish node, but only know +% it is more likely to be on that off +ev2 = ev; +infeng = enter_evidence(infeng, ev2, 'soft_evidence_nodes', [F T], 'soft_evidence', {[0.3 0.7]'}); +end + + +learnednet = learn_params_dbn_em(infeng, {evidence}, 'max_iter', 5); + +disp('real model') +disp_map_hhmm(realnet) + +disp('learned model') +disp_map_hhmm(learnednet) + +disp('rnd model') +disp_map_hhmm(rndnet) + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_map_hhmm.m new file mode 100644 index 00000000..7b077ddb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_map_hhmm.m @@ -0,0 +1,181 @@ +function bnet = mk_map_hhmm(varargin) + +% p is the prob of a successful move (defines the reliability of motors) +p = 1; +obs_model = 'unique'; + +for i=1:2:length(varargin) + switch varargin{i}, + case 'p', p = varargin{i+1}; + case 'obs_model', obs_model = varargin{i+1}; + end +end + + +q = 1-p; +unique_obs = strcmp(obs_model, 'unique'); + +% assign numbers to the nodes in topological order +U = 1; A = 2; C = 3; F = 4; +if unique_obs + onodes = 5; +else + N = 5; E = 6; S = 7; W = 8; % north, east, south, west + onodes = [N E S W]; +end + +% create graph structure + +ss = 4 + length(onodes); % slice size +intra = zeros(ss,ss); +intra(U,F)=1; +intra(A,[C F onodes])=1; +intra(C,[F onodes])=1; + +inter = zeros(ss,ss); +inter(U,[A C])=1; +inter(A,[A C])=1; +inter(F,[A C])=1; +inter(C,C)=1; + +% node sizes +ns = zeros(1,ss); +ns(U) = 2; % left/right +ns(A) = 2; +ns(C) = 3; +ns(F) = 2; +if unique_obs + ns(onodes) = 5; % we will assign each state a unique symbol +else + ns(onodes) = 2; +end +l = 1; r = 2; % left/right +L = 1; R = 2; + +% Make the DBN +bnet = mk_dbn(intra, inter, ns, 'observed', onodes); +eclass = bnet.equiv_class; + + + +% Define CPDs for slice 1 +% We clamp all the CPDs that are not tied, +% since we cannot learn them from a single sequence. + +% uniform probs over actions (the input could be chosen from a policy) +bnet.CPD{eclass(U,1)} = tabular_CPD(bnet, U, 'CPT', mk_stochastic(ones(ns(U),1)), ... + 'adjustable', 0); + +% uniform probs over starting abstract state +bnet.CPD{eclass(A,1)} = tabular_CPD(bnet, A, 'CPT', mk_stochastic(ones(ns(A),1)), ... + 'adjustable', 0); + +% Uniform probs over starting concrete state, modulo the fact +% that corridor 2 is only of length 2. +CPT = zeros(ns(A), ns(C)); % CPT(i,j) = P(C starts in j | A=i) +CPT(1, :) = [1/3 1/3 1/3]; +CPT(2, :) = [1/2 1/2 0]; +bnet.CPD{eclass(C,1)} = tabular_CPD(bnet, C, 'CPT', CPT, 'adjustable', 0); + +% Termination probs +CPT = zeros(ns(U), ns(A), ns(C), ns(F)); +CPT(r,1,1,:) = [1 0]; +CPT(r,1,2,:) = [1 0]; +CPT(r,1,3,:) = [q p]; +CPT(r,2,1,:) = [1 0]; +CPT(r,2,2,:) = [q p]; +CPT(l,1,1,:) = [q p]; +CPT(l,1,2,:) = [1 0]; +CPT(l,1,3,:) = [1 0]; +CPT(l,2,1,:) = [q p]; +CPT(l,2,2,:) = [1 0]; + +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT); + + +% Observation model +if unique_obs + CPT = zeros(ns(A), ns(C), 5); + CPT(1,1,1)=1; % Theo state 4 + CPT(1,2,2)=1; % Theo state 5 + CPT(1,3,3)=1; % Theo state 6 + CPT(2,1,4)=1; % Theo state 9 + CPT(2,2,5)=1; % Theo state 10 + %CPT(2,3,:) undefined + O = onodes(1); + bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT); +else + % north/east/south/west can see wall (1) or opening (2) + CPT = zeros(ns(A), ns(C), 2); + CPT(:,:,1) = q; + CPT(:,:,2) = p; + bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', CPT); + bnet.CPD{eclass(E,1)} = tabular_CPD(bnet, E, 'CPT', CPT); + CPT = zeros(ns(A), ns(C), 2); + CPT(:,:,1) = p; + CPT(:,:,2) = q; + bnet.CPD{eclass(S,1)} = tabular_CPD(bnet, S, 'CPT', CPT); + bnet.CPD{eclass(N,1)} = tabular_CPD(bnet, N, 'CPT', CPT); +end + +% Define the CPDs for slice 2 + +% Abstract + +% Since the top level never resets, the starting distribution is irrelevant: +% A2 will be determined by sampling from transmat(A1,:). +% But the code requires we specify it anyway; we make it all 0s, a dummy value. +startprob = zeros(ns(U), ns(A)); + +transmat = zeros(ns(U), ns(A), ns(A)); +transmat(R,1,:) = [q p]; +transmat(R,2,:) = [0 1]; +transmat(L,1,:) = [1 0]; +transmat(L,2,:) = [p q]; + +% Qps are the parents we condition the parameters on, in this case just +% the past action. +bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ... + 'startprob', startprob, 'transprob', transmat); + + + +% Concrete + +transmat = zeros(ns(C), ns(U), ns(A), ns(C)); +transmat(1,r,1,:) = [q p 0.0]; +transmat(2,r,1,:) = [0.0 q p]; +transmat(3,r,1,:) = [0.0 0.0 1.0]; +transmat(1,r,2,:) = [q p 0.0]; +transmat(2,r,2,:) = [0.0 1.0 0.0]; +% +transmat(1,l,1,:) = [1.0 0.0 0.0]; +transmat(2,l,1,:) = [p q 0.0]; +transmat(3,l,1,:) = [0.0 p q]; +transmat(1,l,2,:) = [1.0 0.0 0.0]; +transmat(2,l,2,:) = [p q 0.0]; + +% Add a new dimension for A(t-1), by copying old vals, +% so the matrix is the same size as startprob + + +transmat = reshape(transmat, [ns(C) ns(U) ns(A) 1 ns(C)]); +transmat = repmat(transmat, [1 1 1 ns(A) 1]); + +% startprob(C(t-1), U(t-1), A(t-1), A(t), C(t)) +startprob = zeros(ns(C), ns(U), ns(A), ns(A), ns(C)); +startprob(1,L,1,1,:) = [1.0 0.0 0.0]; +startprob(3,R,1,2,:) = [1.0 0.0 0.0]; +startprob(3,R,1,1,:) = [0.0 0.0 1.0]; +% +startprob(1,L,2,1,:) = [0.0 0.0 010]; +startprob(2,L,2,1,:) = [1.0 0.0 0.0]; +startprob(2,R,2,2,:) = [0.0 1.0 0.0]; + +% want transmat(U,A,C,At,Ct), ie. in topo order +transmat = permute(transmat, [2 3 1 4 5]); +startprob = permute(startprob, [2 3 1 4 5]); +bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ... + 'startprob', startprob, 'transprob', transmat); + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_rnd_map_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_rnd_map_hhmm.m new file mode 100644 index 00000000..76b06fc7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/mk_rnd_map_hhmm.m @@ -0,0 +1,73 @@ +function bnet = mk_rnd_map_hhmm(varargin) + +% We copy the deterministic structure of the real HHMM, +% but randomize the probabilities of the adjustable CPDs. +% The key trick is that 0s in the real HHMM remain 0 +% even when multiplied by a randon number. + +obs_model = 'unique'; + +for i=1:2:length(varargin) + switch varargin{i}, + case 'obs_model', obs_model = varargin{i+1}; + end +end + + +unique_obs = strcmp(obs_model, 'unique'); + +psuccess = 0.9; +% must be less than 1, so that pfail > 0 +% otherwise we copy too many 0s +bnet = mk_map_hhmm('p', psuccess, 'obs_model', obs_model); +ns = bnet.node_sizes; +ss = bnet.nnodes_per_slice; + +U = 1; A = 2; C = 3; F = 4; +%unique_obs = (bnet.nnodes_per_slice == 5); +if unique_obs + onodes = 5; +else + north = 5; east = 6; south = 7; west = 8; + onodes = [north east south west]; +end + +eclass = bnet.equiv_class; +S=struct(bnet.CPD{eclass(F,1)}); +CPT = mk_stochastic(rand(size(S.CPT)) .* S.CPT); +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT); + + +% Observation model +if unique_obs + CPT = zeros(ns(A), ns(C), 5); + CPT(1,1,1)=1; % Theo state 4 + CPT(1,2,2)=1; % Theo state 5 + CPT(1,3,3)=1; % Theo state 6 + CPT(2,1,4)=1; % Theo state 9 + CPT(2,2,5)=1; % Theo state 10 + %CPT(2,3,:) undefined + O = onodes(1); + bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT); +else + for i=[north east south west] + CPT = mk_stochastic(rand(ns(A), ns(C), 2)); + bnet.CPD{eclass(i,1)} = tabular_CPD(bnet, i, 'CPT', CPT); + end +end + +% Define the CPDs for slice 2 + +startprob = zeros(ns(U), ns(A)); +S = struct(bnet.CPD{eclass(A,2)}); +transprob = mk_stochastic(rand(size(S.transprob)) .* S.transprob); +bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ... + 'startprob', startprob, 'transprob', transprob); + +S = struct(bnet.CPD{eclass(C,2)}); +transprob = mk_stochastic(rand(size(S.transprob)) .* S.transprob); +startprob = mk_stochastic(rand(size(S.startprob)) .* S.startprob); +bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ... + 'startprob', startprob, 'transprob', transprob); + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/sample_from_map.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/sample_from_map.m new file mode 100644 index 00000000..816b741e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Map/sample_from_map.m @@ -0,0 +1,41 @@ +if 0 +% Generate some sample paths + +bnet = mk_map_hhmm('p', 1); +% assign numbers to the nodes in topological order +U = 1; A = 2; C = 3; F = 4; O = 5; + + +seed = 0; +rand('state', seed); +randn('state', seed); + +% control policy = sweep right then left +T = 10; +ss = 5; +ev = cell(ss, T); +ev(U,:) = num2cell([R*ones(1,5) L*ones(1,5)]); + +% fix initial conditions to be in left most state +ev{A,1} = 1; +ev{C,1} = 1; +evidence = sample_dbn(bnet, 'length', T, 'evidence', ev) + + +% Now do same but with noisy actuators + +bnet = mk_map_hhmm('p', 0.8); +evidence = sample_dbn(bnet, 'length', T, 'evidence', ev) + +end + +% Now do same but with 4 observations per slice + +bnet = mk_map_hhmm('p', 0.8, 'obs_model', 'four'); +ss = bnet.nnodes_per_slice; + +ev = cell(ss, T); +ev(U,:) = num2cell([R*ones(1,5) L*ones(1,5)]); +ev{A,1} = 1; +ev{C,1} = 1; +evidence = sample_dbn(bnet, 'length', T, 'evidence', ev) diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Entries new file mode 100644 index 00000000..c4379581 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Entries @@ -0,0 +1,6 @@ +/letter2num.m/1.1.1.1/Fri Nov 22 23:10:20 2002// +/mgram1.m/1.1.1.1/Fri Nov 22 23:59:00 2002// +/mgram2.m/1.1.1.1/Tue Nov 26 22:04:24 2002// +/mgram3.m/1.1.1.1/Tue Nov 26 22:14:10 2002// +/num2letter.m/1.1.1.1/Fri Nov 22 23:07:40 2002// +D/Old//// diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Repository new file mode 100644 index 00000000..5ede715e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Mgram diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Entries new file mode 100644 index 00000000..07b688e2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Entries @@ -0,0 +1,2 @@ +/mgram2.m/1.1.1.1/Sat Nov 23 00:44:34 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Repository new file mode 100644 index 00000000..ba3ae6df --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Mgram/Old diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/mgram2.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/mgram2.m new file mode 100644 index 00000000..719a3167 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/Old/mgram2.m @@ -0,0 +1,191 @@ +% like mgram1, except we use a durational HMM instead of an HHMM2 + +past = 0; + +words = {'the', 't', 'h', 'e'}; +data = 'the'; +nwords = length(words); +word_len = zeros(1, nwords); +word_prob = normalise(ones(1,nwords)); +word_logprob = log(word_prob); +for wi=1:nwords + word_len(wi)=length(words{wi}); +end +D = max(word_len); + + +alphasize = 26*2; +data = letter2num(data); +T = length(data); + +% node numbers +W = 1; % top level state = word id +L = 2; % bottom level state = letter position within word +F = 3; +O = 4; + +ss = 4; +intra = zeros(ss,ss); +intra(W,[F L O])=1; +intra(L,[O F])=1; + +inter = zeros(ss,ss); +inter(W,W)=1; +inter(L,L)=1; +inter(F,[W L O])=1; + +% node sizes +ns = zeros(1,ss); +ns(W) = nwords; +ns(L) = D; +ns(F) = 2; +ns(O) = alphasize; +ns2 = [ns ns]; + +% Make the DBN +bnet = mk_dbn(intra, inter, ns, 'observed', O); +eclass = bnet.equiv_class; + +% uniform start distrib over words, uniform trans mat +Wstart = normalise(ones(1,nwords)); +Wtrans = mk_stochastic(ones(nwords,nwords)); + +% always start in state d = length(word) for each bottom level HMM +Lstart = zeros(nwords, D); +for i=1:nwords + l = length(words{i}); + Lstart(i,l)=1; +end + +% make downcounters +RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob +Ltrans = repmat(RLtrans, [1 1 nwords]); + +% Finish when downcoutner = 1 +Fprob = zeros(nwords, D, 2); +Fprob(:,1,2)=1; +Fprob(:,2:end,1)=1; + + +% Define CPDs for slice +bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart); +bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart); +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob); + + +% Define CPDs for slice 2 +bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart, 'transprob', Wtrans); +bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans); + + +if 0 +% To test it is generating correctly, we create an artificial +% observation process that capitalizes at the start of a new segment +% Oprob(Ft-1,Qt,Dt,Yt) +Oprob = zeros(2,nwords,D,alphasize); +Oprob(1,1,3,letter2num('t'),1)=1; +Oprob(1,1,2,letter2num('h'),1)=1; +Oprob(1,1,1,letter2num('e'),1)=1; +Oprob(2,1,3,letter2num('T'),1)=1; +Oprob(2,1,2,letter2num('H'),1)=1; +Oprob(2,1,1,letter2num('E'),1)=1; +Oprob(1,2,1,letter2num('a'),1)=1; +Oprob(2,2,1,letter2num('A'),1)=1; +Oprob(1,3,1,letter2num('b'),1)=1; +Oprob(2,3,1,letter2num('B'),1)=1; +Oprob(1,4,1,letter2num('c'),1)=1; +Oprob(2,4,1,letter2num('C'),1)=1; + +% Oprob1(Qt,Dt,Yt) +Oprob1 = zeros(nwords,D,alphasize); +Oprob1(1,3,letter2num('t'),1)=1; +Oprob1(1,2,letter2num('h'),1)=1; +Oprob1(1,1,letter2num('e'),1)=1; +Oprob1(2,1,letter2num('a'),1)=1; +Oprob1(3,1,letter2num('b'),1)=1; +Oprob1(4,1,letter2num('c'),1)=1; + +bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob); +bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1); + +evidence = cell(ss,T); +%evidence{W,1}=1; +sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence)); +str = num2letter(sample(4,:)) +end + + + + +[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob); +% obslik(j,t,d) +softCPDpot = cell(ss,T); +ens = ns; +ens(O)=1; +ens2 = [ens ens]; +for t=2:T + dom = [F W+ss L+ss O+ss]; + % tab(Ft-1, Q2, Dt) + tab = ones(2, nwords, D); + if past + tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1 + tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment + else + for d=1:max(1,min(D,T+1-t)) + tab(2,:,d) = squeeze(obslik(:,t+d-1,d)); + end + end + softCPDpot{O,t} = dpot(dom, ens2(dom), tab); +end +t = 1; +dom = [W L O]; +% tab(Q2, Dt) +tab = ones(nwords, D); +if past + tab = squeeze(obslik(:,t,:)); +else + for d=1:min(D,T-t) + tab(:,d) = squeeze(obslik(:,t+d-1,d)); + end +end +softCPDpot{O,t} = dpot(dom, ens(dom), tab); + + +%bnet.observed = []; +% uniformative observations +%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize))); +%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize))); + +engine = jtree_dbn_inf_engine(bnet); +evidence = cell(ss,T); +% we add dummy data to O to force its effective size to be 1. +% The actual values have already been incorporated into softCPDpot +evidence(O,:) = num2cell(ones(1,T)); +[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot); + + +%evidence(F,:) = num2cell(2*ones(1,T)); +%[engine, ll_dbn] = enter_evidence(engine, evidence); + + +gamma = zeros(nwords, T); +for t=1:T + m = marginal_nodes(engine, [W F], t); + gamma(:,t) = m.T(:,2); +end + +gamma + +xidbn = zeros(nwords, nwords); +for t=1:T-1 + m = marginal_nodes(engine, [W F W+ss], t); + xidbn = xidbn + squeeze(m.T(:,2,:)); +end + +% thee +% xidbn(1,4) = 0.9412 the->e +% (2,3)=0.0588 t->h +% (3,4)=0.0588 h-e +% (4,4)=0.0588 e-e + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/letter2num.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/letter2num.m new file mode 100644 index 00000000..f4e3f1d9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/letter2num.m @@ -0,0 +1,12 @@ +function n = letter2num(l) + +% map a-z to 1:26 and A-Z to 27:52 +punct_code = [32:47 58:64 91:96 123:126]; +digits_code = 48:57; +upper_code = 65:90; +lower_code = 97:122; + +c = double(l); +n = c-96; +ndx = find(n <= 0); % upper case +n(ndx) = c(ndx) - 64 + 26; diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram1.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram1.m new file mode 100644 index 00000000..52ca472e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram1.m @@ -0,0 +1,116 @@ +% a multigram is a degenerate 2HHMM where the bottom level HMMs emit deterministic strings +% and the the top level abstract states are independent of each other +% cf. HSMM/test_mgram2 + +words = {'the', 't', 'h', 'e'}; +data = 'the'; +nwords = length(words); +word_len = zeros(1, nwords); +word_prob = normalise(ones(1,nwords)); +word_logprob = log(word_prob); +for wi=1:nwords + word_len(wi)=length(words{wi}); +end +D = max(word_len); + +alphasize = 26; +data = letter2num(data); +T = length(data); + +% node numbers +W = 1; % top level state = word id +L = 2; % bottom level state = letter position within word +F = 3; +O = 4; + +ss = 4; +intra = zeros(ss,ss); +intra(W,[F L O])=1; +intra(L,[O F])=1; + +inter = zeros(ss,ss); +inter(W,W)=1; +inter(L,L)=1; +inter(F,[W L])=1; + +% node sizes +ns = zeros(1,ss); +ns(W) = nwords; +ns(L) = D; +ns(F) = 2; +ns(O) = alphasize; + + +% Make the DBN +bnet = mk_dbn(intra, inter, ns, 'observed', O); +eclass = bnet.equiv_class; + + + +% uniform start distrib over words, uniform trans mat +Wstart = normalise(ones(1,nwords)); +Wtrans = mk_stochastic(ones(nwords,nwords)); + +% always start in state 1 for each bottom level HMM +delta1_start = zeros(1, D); +delta1_start(1) = 1; +Lstart = repmat(delta1_start, nwords, 1); +LRtrans = mk_leftright_transmat(D, 0); % 0 self loop prob +Ltrans = repmat(LRtrans, [1 1 nwords]); + +% Finish in the last letter of each word +Fprob = zeros(nwords, D, 2); +Fprob(:,:,1)=1; +for i=1:nwords + Fprob(i,length(words{i}),2)=1; + Fprob(i,length(words{i}),1)=0; +end + +% Each state uniquely emits a letter +Oprob = zeros(nwords, D, alphasize); +for i=1:nwords + for l=1:length(words{i}) + a = double(words{i}(l))-96; + Oprob(i,l,a)=1; + end +end + + +% Define CPDs for slice +bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart); +bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart); +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob); +bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob); + +% Define CPDs for slice 2 +bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart, 'transprob', Wtrans); +bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans); + +evidence = cell(ss,T); +evidence{W,1}=1; +sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence)); +str = lower(sample(4,:)) + +engine = jtree_dbn_inf_engine(bnet); +evidence = cell(ss,T); +evidence(O,:) = num2cell(data); +[engine, ll_dbn] = enter_evidence(engine, evidence); + +gamma = zeros(nwords, T); +for t=1:T + m = marginal_nodes(engine, [W F], t); + gamma(:,t) = m.T(:,2); +end +gamma + +xidbn = zeros(nwords, nwords); +for t=1:T-1 + m = marginal_nodes(engine, [W F W+ss], t); + xidbn = xidbn + squeeze(m.T(:,2,:)); +end + +% thee +% xidbn(1,4) = 0.9412 the->e +% (2,3)=0.0588 t->h +% (3,4)=0.0588 h-e +% (4,4)=0.0588 e-e diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram2.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram2.m new file mode 100644 index 00000000..c61f855a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram2.m @@ -0,0 +1,200 @@ +% Like a durational HMM, except we use soft evidence on the observed nodes. +% Should give the same results as HSMM/test_mgram2. + +past = 1; +% If past=1, P(Yt|Qt=j,Dt=d) = P(y_{t-d+1:t}|j) +% If past=0, P(Yt|Qt=j,Dt=d) = P(y_{t:t+d-1}|j) - future evidence + +words = {'the', 't', 'h', 'e'}; +data = 'the'; +nwords = length(words); +word_len = zeros(1, nwords); +word_prob = normalise(ones(1,nwords)); +word_logprob = log(word_prob); +for wi=1:nwords + word_len(wi)=length(words{wi}); +end +D = max(word_len); + + +alphasize = 26*2; +data = letter2num(data); +T = length(data); + +% node numbers +W = 1; % top level state = word id +L = 2; % bottom level state = letter position within word +F = 3; +O = 4; + +ss = 4; +intra = zeros(ss,ss); +intra(W,[F L O])=1; +intra(L,[O F])=1; + +inter = zeros(ss,ss); +inter(W,W)=1; +inter(L,L)=1; +inter(F,[W L O])=1; + +% node sizes +ns = zeros(1,ss); +ns(W) = nwords; +ns(L) = D; +ns(F) = 2; +ns(O) = alphasize; +ns2 = [ns ns]; + +% Make the DBN +bnet = mk_dbn(intra, inter, ns, 'observed', O); +eclass = bnet.equiv_class; + +% uniform start distrib over words, uniform trans mat +Wstart = normalise(ones(1,nwords)); +Wtrans = mk_stochastic(ones(nwords,nwords)); +%Wtrans = ones(nwords,nwords); + +% always start in state d = length(word) for each bottom level HMM +Lstart = zeros(nwords, D); +for i=1:nwords + l = length(words{i}); + Lstart(i,l)=1; +end + +% make downcounters +RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob +Ltrans = repmat(RLtrans, [1 1 nwords]); + +% Finish when downcoutner = 1 +Fprob = zeros(nwords, D, 2); +Fprob(:,1,2)=1; +Fprob(:,2:end,1)=1; + + +% Define CPDs for slice 1 +bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart); +bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart); +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob); + + +% Define CPDs for slice 2 +bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart, 'transprob', Wtrans); +bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans); + + +if 0 +% To test it is generating correctly, we create an artificial +% observation process that capitalizes at the start of a new segment +% Oprob(Ft-1,Qt,Dt,Yt) +Oprob = zeros(2,nwords,D,alphasize); +Oprob(1,1,3,letter2num('t'),1)=1; +Oprob(1,1,2,letter2num('h'),1)=1; +Oprob(1,1,1,letter2num('e'),1)=1; +Oprob(2,1,3,letter2num('T'),1)=1; +Oprob(2,1,2,letter2num('H'),1)=1; +Oprob(2,1,1,letter2num('E'),1)=1; +Oprob(1,2,1,letter2num('a'),1)=1; +Oprob(2,2,1,letter2num('A'),1)=1; +Oprob(1,3,1,letter2num('b'),1)=1; +Oprob(2,3,1,letter2num('B'),1)=1; +Oprob(1,4,1,letter2num('c'),1)=1; +Oprob(2,4,1,letter2num('C'),1)=1; + +% Oprob1(Qt,Dt,Yt) +Oprob1 = zeros(nwords,D,alphasize); +Oprob1(1,3,letter2num('t'),1)=1; +Oprob1(1,2,letter2num('h'),1)=1; +Oprob1(1,1,letter2num('e'),1)=1; +Oprob1(2,1,letter2num('a'),1)=1; +Oprob1(3,1,letter2num('b'),1)=1; +Oprob1(4,1,letter2num('c'),1)=1; + +bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob); +bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1); + +evidence = cell(ss,T); +%evidence{W,1}=1; +sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence)); +str = num2letter(sample(4,:)) +end + + +if 1 + +[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob); +% obslik(j,t,d) +softCPDpot = cell(ss,T); +ens = ns; +ens(O)=1; +ens2 = [ens ens]; +for t=2:T + dom = [F W+ss L+ss O+ss]; + % tab(Ft-1, Q2, Dt) + tab = ones(2, nwords, D); + if past + tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1 + %tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment + for d=1:min(t,D) + tab(2,:,d) = squeeze(obslik(:,t,d)); + end + else + for d=1:max(1,min(D,T+1-t)) + tab(2,:,d) = squeeze(obslik(:,t+d-1,d)); + end + end + softCPDpot{O,t} = dpot(dom, ens2(dom), tab); +end +t = 1; +dom = [W L O]; +% tab(Q2, Dt) +tab = ones(nwords, D); +if past + %tab = squeeze(obslik(:,t,:)); + tab(:,1) = squeeze(obslik(:,t,1)); +else + for d=1:min(D,T-t) + tab(:,d) = squeeze(obslik(:,t+d-1,d)); + end +end +softCPDpot{O,t} = dpot(dom, ens(dom), tab); + + +%bnet.observed = []; +% uniformative observations +%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize))); +%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize))); + +engine = jtree_dbn_inf_engine(bnet); +evidence = cell(ss,T); +% we add dummy data to O to force its effective size to be 1. +% The actual values have already been incorporated into softCPDpot +evidence(O,:) = num2cell(ones(1,T)); +[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot); + + +%evidence(F,:) = num2cell(2*ones(1,T)); +%[engine, ll_dbn] = enter_evidence(engine, evidence); + + +gamma = zeros(nwords, T); +for t=1:T + m = marginal_nodes(engine, [W F], t); + gamma(:,t) = m.T(:,2); +end + +gamma + +xidbn = zeros(nwords, nwords); +for t=1:T-1 + m = marginal_nodes(engine, [W F W+ss], t); + xidbn = xidbn + squeeze(m.T(:,2,:)); +end + +% thee +% xidbn(1,4) = 0.9412 the->e +% (2,3)=0.0588 t->h +% (3,4)=0.0588 h-e +% (4,4)=0.0588 e-e + + +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram3.m new file mode 100644 index 00000000..49228584 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/mgram3.m @@ -0,0 +1,235 @@ +% like mgram2, except we unroll the DBN so we can use smaller +% state spaces for the early duration nodes: +% the state spaces are D1 in {1}, D2 in {1,2} + +past = 1; + +words = {'the', 't', 'h', 'e'}; +data = 'the'; +nwords = length(words); +word_len = zeros(1, nwords); +word_prob = normalise(ones(1,nwords)); +word_logprob = log(word_prob); +for wi=1:nwords + word_len(wi)=length(words{wi}); +end +D = max(word_len); + + +alphasize = 26*2; +data = letter2num(data); +T = length(data); + +% node numbers +W = 1; % top level state = word id +L = 2; % bottom level state = letter position within word +F = 3; +O = 4; + +ss = 4; +intra = zeros(ss,ss); +intra(W,[F L O])=1; +intra(L,[O F])=1; + +inter = zeros(ss,ss); +inter(W,W)=1; +inter(L,L)=1; +inter(F,[W L O])=1; + +T = 3; +dag = unroll_dbn_topology(intra, inter, T); + +% node sizes +ns = zeros(1,ss); +ns(W) = nwords; +ns(L) = D; +ns(F) = 2; +ns(O) = alphasize; +ns = repmat(ns(:), [1 T]); +for d=1:D + ns(d,L)=d; % max duration +end +ns = ns(:); + +% Equiv class in brackets for D=3 +% The Lt's are not tied until t>=D, since they have different sizes. +% W1 and W2 are not tied since they have different parent sets. + +% W1 (1) W2 (5) W3 (5) W4 (5) +% L1 (2) L2 (6) L3 (7) L4 (7) +% F1 (3) F2 (3) F3 (4) F3 (4) +% O1 (4) O2 (4) O2 (4) O4 (4) + +% Since we are not learning, we can dispense with tying + +% Make the bnet +Wnodes = unroll_set(W, ss, T); +Lnodes = unroll_set(L, ss, T); +Fnodes = unroll_set(F, ss, T); +Onodes = unroll_set(O, ss, T); + +bnet = mk_bnet(dag, ns); +eclass = bnet.equiv_class; + +% uniform start distrib over words, uniform trans mat +Wstart = normalise(ones(1,nwords)); +Wtrans = mk_stochastic(ones(nwords,nwords)); +bnet.CPD{eclass(Wnodes(1))} = tabular_CPD(bnet, Wnodes(1), 'CPT', Wstart); +for t=2:T +bnet.CPD{eclass(Wnodes(t))} = hhmmQ_CPD(bnet, Wnodes(t), 'Fbelow', Fnodes(t-1), ... + 'startprob', Wstart, 'transprob', Wtrans); +end + +% always start in state d = length(word) for each bottom level HMM +% and then count down +% make downcounters +RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob +Ltrans = repmat(RLtrans, [1 1 nwords]); + +for t=1:T + Lstart = zeros(nwords, min(t,D)); + for i=1:nwords + l = length(words{i}); + Lstart(i,l)=1; + if d==1 + bnet.CPD{eclass(Lnodes(1))} = tabular_CPD(bnet, Lnodes(1), 'CPT', Lstart); + else + bnet.CPD{eclass(Lnodes(t))} = hhmmQ_CPD(bnet, Lnodes(t), 'Fself', Fnodes(t-1), 'Qps', Wnodes(t), ... + 'startprob', Lstart, 'transprob', Ltrans); + end + end +end + + +% Finish when downcoutner = 1 +Fprob = zeros(nwords, D, 2); +Fprob(:,1,2)=1; +Fprob(:,2:end,1)=1; + + +% Define CPDs for slice +bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart); +bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart); +bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob); + + +% Define CPDs for slice 2 +bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart, 'transprob', Wtrans); +bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans); + + +if 0 +% To test it is generating correctly, we create an artificial +% observation process that capitalizes at the start of a new segment +% Oprob(Ft-1,Qt,Dt,Yt) +Oprob = zeros(2,nwords,D,alphasize); +Oprob(1,1,3,letter2num('t'),1)=1; +Oprob(1,1,2,letter2num('h'),1)=1; +Oprob(1,1,1,letter2num('e'),1)=1; +Oprob(2,1,3,letter2num('T'),1)=1; +Oprob(2,1,2,letter2num('H'),1)=1; +Oprob(2,1,1,letter2num('E'),1)=1; +Oprob(1,2,1,letter2num('a'),1)=1; +Oprob(2,2,1,letter2num('A'),1)=1; +Oprob(1,3,1,letter2num('b'),1)=1; +Oprob(2,3,1,letter2num('B'),1)=1; +Oprob(1,4,1,letter2num('c'),1)=1; +Oprob(2,4,1,letter2num('C'),1)=1; + +% Oprob1(Qt,Dt,Yt) +Oprob1 = zeros(nwords,D,alphasize); +Oprob1(1,3,letter2num('t'),1)=1; +Oprob1(1,2,letter2num('h'),1)=1; +Oprob1(1,1,letter2num('e'),1)=1; +Oprob1(2,1,letter2num('a'),1)=1; +Oprob1(3,1,letter2num('b'),1)=1; +Oprob1(4,1,letter2num('c'),1)=1; + +bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob); +bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1); + +evidence = cell(ss,T); +%evidence{W,1}=1; +sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence)); +str = num2letter(sample(4,:)) +end + + + + +[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob); +% obslik(j,t,d) +softCPDpot = cell(ss,T); +ens = ns; +ens(O)=1; +ens2 = [ens ens]; +for t=2:T + dom = [F W+ss L+ss O+ss]; + % tab(Ft-1, Q2, Dt) + tab = ones(2, nwords, D); + if past + tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1 + %tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment + for d=1:min(t,D) + tab(2,:,d) = squeeze(obslik(:,t,d)); + end + else + for d=1:max(1,min(D,T+1-t)) + tab(2,:,d) = squeeze(obslik(:,t+d-1,d)); + end + end + softCPDpot{O,t} = dpot(dom, ens2(dom), tab); +end +t = 1; +dom = [W L O]; +% tab(Q2, Dt) +tab = ones(nwords, D); +if past + %tab = squeeze(obslik(:,t,:)); + tab(:,1) = squeeze(obslik(:,t,1)); +else + for d=1:min(D,T-t) + tab(:,d) = squeeze(obslik(:,t+d-1,d)); + end +end +softCPDpot{O,t} = dpot(dom, ens(dom), tab); + + +%bnet.observed = []; +% uniformative observations +%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize))); +%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize))); + +engine = jtree_dbn_inf_engine(bnet); +evidence = cell(ss,T); +% we add dummy data to O to force its effective size to be 1. +% The actual values have already been incorporated into softCPDpot +evidence(O,:) = num2cell(ones(1,T)); +[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot); + + +%evidence(F,:) = num2cell(2*ones(1,T)); +%[engine, ll_dbn] = enter_evidence(engine, evidence); + + +gamma = zeros(nwords, T); +for t=1:T + m = marginal_nodes(engine, [W F], t); + gamma(:,t) = m.T(:,2); +end + +gamma + +xidbn = zeros(nwords, nwords); +for t=1:T-1 + m = marginal_nodes(engine, [W F W+ss], t); + xidbn = xidbn + squeeze(m.T(:,2,:)); +end + +% thee +% xidbn(1,4) = 0.9412 the->e +% (2,3)=0.0588 t->h +% (3,4)=0.0588 h-e +% (4,4)=0.0588 e-e + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/num2letter.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/num2letter.m new file mode 100644 index 00000000..139b81fa --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Mgram/num2letter.m @@ -0,0 +1,10 @@ +function l = num2letter(n) + +% map 1:26 to a-z and 27:52 to A-Z +punct_code = [32:47 58:64 91:96 123:126]; +digits_code = 48:57; +upper_code = 65:90; +lower_code = 97:122; + +letters = [char(lower_code) char(upper_code)]; +l = letters(n); diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries new file mode 100644 index 00000000..93258c2a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries @@ -0,0 +1,5 @@ +/fixed_args_mk_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/learn_motif_hhmm.m/1.1.1.1/Tue Jul 2 22:56:14 2002// +/mk_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository new file mode 100644 index 00000000..5062cfd4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Motif diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m new file mode 100644 index 00000000..ce4e288b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m @@ -0,0 +1,99 @@ +function bnet = fixed_args_mk_motif_hhmm(motif_length, motif_pattern, background_char) +% +% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH) +% Make the following HHMM +% +% S2 <----------------------> S1 +% | | +% | | +% M1 -> M2 -> M3 -> end B1 -> end +% +% where Mi represents the i'th letter in the motif +% and B is the background state. +% Si chooses between running the motif or the background. +% The Si and B states have self loops (not shown). +% +% The transition params are defined to respect the above topology. +% The background is uniform; each motif state has a random obs. distribution. +% +% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN) +% In this case, we make the motif submodel deterministically +% emit the motif pattern. +% +% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN, BACKGROUND_CHAR) +% In this case, we make the background submodel +% deterministically emit the specified character (to make the pattern +% easier to see). + +if nargin < 2, motif_pattern = []; end +if nargin < 3, background_char = []; end + +chars = ['a', 'c', 'g', 't']; +Osize = length(chars); + +motif_length = length(motif_pattern); +Qsize = [2 motif_length]; +Qnodes = 1:2; +D = 2; +transprob = cell(1,D); +termprob = cell(1,D); +startprob = cell(1,D); + +% startprob{d}(k,j), startprob{1}(1,j) +% transprob{d}(i,k,j), transprob{1}(i,j) +% termprob{d}(k,j) + + +% LEVEL 1 + +startprob{1} = zeros(1, 2); +startprob{1} = [1 0]; % always start in the background model + +% When in the background state, we stay there with high prob +% When in the motif state, we immediately return to the background state. +transprob{1} = [0.8 0.2; + 1.0 0.0]; + + +% LEVEL 2 +startprob{2} = 'leftstart'; % both submodels start in substate 1 +transprob{2} = zeros(motif_length, 2, motif_length); +termprob{2} = zeros(2, motif_length); + +% In the background model, we only use state 1. +transprob{2}(1,1,1) = 1; % self loop +termprob{2}(1,1) = 0.2; % prob transition to end state + +% Motif model +transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); % no self loops +termprob{2}(2,end) = 1.0; % last state immediately terminates + + +% OBS LEVEl + +obsprob = zeros([Qsize Osize]); +if isempty(background_char) + % uniform background model + obsprob(1,1,:) = normalise(ones(Osize,1)); +else + % deterministic background model (easy to see!) + m = find(chars==background_char); + obsprob(1,1,m) = 1.0; +end + +if gen_motif + % initialise with true motif (cheating) + for i=1:motif_length + m = find(chars == motif_pattern(i)); + obsprob(2,i,m) = 1.0; + end +else + obsprob(2,:,:) = mk_stochastic(ones(motif_length, Osize)); +end + +Oargs = {'CPT', obsprob}; + +[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ... + 'Oargs', Oargs, 'Ops', Qnodes(1:2), ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m new file mode 100644 index 00000000..54553460 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m @@ -0,0 +1,75 @@ + +seed = 0; +rand('state', seed); +randn('state', seed); + +chars = ['a', 'c', 'g', 't']; +motif = 'accca'; +motif_length = length(motif); +motif_code = zeros(1, motif_length); +for i=1:motif_length + motif_code(i) = find(chars == motif(i)); +end + +[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_length', length(motif)); +%[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_pattern', motif); +ss = bnet_init.nnodes_per_slice; + + + +% We generate a training set by creating uniform sequences, +% and inserting a single motif at a random location. +ntrain = 100; +T = 20; +cases = cell(1, ntrain); + +if 1 + % uniform background + background_dist = normalise(ones(1, length(chars))); +end +if 0 + % use a constant background + background_dist = zeros(1, length(chars)); + m = find(chars=='t'); + background_dist(m) = 1.0; +end +if 0 + % use a background skewed away from the motif + p = 0.01; q = (1-(2*p))/2; + background_dist = [p p q q]; +end + +unif_pos = normalise(ones(1, T-length(motif))); +cases = cell(1, ntrain); +data = zeros(1,T); +for i=1:ntrain + data = sample_discrete(background_dist, 1, T); + L = sample_discrete(unif_pos, 1, 1); + data(L:L+length(motif)-1) = motif_code; + cases{i} = cell(ss, T); + cases{i}(Onode,:) = num2cell(data); +end +disp('sample training cases') +for i=1:5 + chars(cell2num(cases{i}(Onode,:))) +end + +engine_init = hmm_inf_engine(bnet_init); + +[bnet_learned, LL, engine_learned] = ... + learn_params_dbn_em(engine_init, cases, 'max_iter', 100, 'thresh', 1e-2); +% 'anneal', 1, 'anneal_rate', 0.7); + +% extract the learned motif profile +eclass = bnet_learned.equiv_class; +CPDO=struct(bnet_learned.CPD{eclass(Onode,1)}); +fprintf('columns = chars, rows = states\n'); +profile_learned = squeeze(CPDO.CPT(2,:,:)) +[m,ndx] = max(profile_learned, [], 2); +map_motif_learned = chars(ndx) +back_learned = squeeze(CPDO.CPT(1,1,:))' +%map_back_learned = chars(argmax(back_learned)) + +CPDO_init = struct(bnet_init.CPD{eclass(Onode,1)}); +profile_init = squeeze(CPDO_init.CPT(2,:,:)); +back_init = squeeze(CPDO_init.CPT(1,1,:))'; diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m new file mode 100644 index 00000000..32980397 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m @@ -0,0 +1,137 @@ +function [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(varargin) +% [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(...) +% +% Make the following HHMM +% +% S2 <----------------------> S1 +% | | +% | | +% M1 -> M2 -> M3 -> end B1 -> end +% +% where Mi represents the i'th letter in the motif +% and B is the background state. +% Si chooses between running the motif or the background. +% The Si and B states have self loops (not shown). +% +% The transition params are defined to respect the above topology. +% The background is uniform; each motif state has a random obs. distribution. +% +% Optional params: +% motif_length - required, unless we specify motif_pattern +% motif_pattern - if specified, we make the motif submodel deterministically +% emit this pattern +% background - if specified, we make the background submodel +% deterministically emit this (makes the motif easier to see!) + + +args = varargin; +nargs = length(args); + +% extract pattern, if any +motif_pattern = []; +for i=1:2:nargs + switch args{i}, + case 'motif_pattern', motif_pattern = args{i+1}; + end +end + +% set defaults +motif_length = length(motif_pattern); +background_char = []; + +% get params +for i=1:2:nargs + switch args{i}, + case 'motif_length', motif_length = args{i+1}; + case 'background', background_char = args{i+1}; + end +end + + +chars = ['a', 'c', 'g', 't']; +Osize = length(chars); + +Qsize = [2 motif_length]; +Qnodes = 1:2; +D = 2; +transprob = cell(1,D); +termprob = cell(1,D); +startprob = cell(1,D); + +% startprob{d}(k,j), startprob{1}(1,j) +% transprob{d}(i,k,j), transprob{1}(i,j) +% termprob{d}(k,j) + + +% LEVEL 1 + +startprob{1} = zeros(1, 2); +startprob{1} = [1 0]; % always start in the background model + +% When in the background state, we stay there with high prob +% When in the motif state, we immediately return to the background state. +transprob{1} = [0.8 0.2; + 1.0 0.0]; + + +% LEVEL 2 +startprob{2} = 'leftstart'; % both submodels start in substate 1 +transprob{2} = zeros(motif_length, 2, motif_length); +termprob{2} = zeros(2, motif_length); + +% In the background model, we only use state 1. +transprob{2}(1,1,1) = 1; % self loop +termprob{2}(1,1) = 0.2; % prob transition to end state + +% Motif model +transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); % no self loops +termprob{2}(2,end) = 1.0; % last state immediately terminates + + +% OBS LEVEl + +obsprob = zeros([Qsize Osize]); +if isempty(background_char) + % uniform background model + %obsprob(1,1,:) = normalise(ones(Osize,1)); + obsprob(1,1,:) = normalise(rand(Osize,1)); +else + % deterministic background model (easy to see!) + m = find(chars==background_char); + obsprob(1,1,m) = 1.0; +end + +if ~isempty(motif_pattern) + % initialise with true motif (cheating) + for i=1:motif_length + m = find(chars == motif_pattern(i)); + obsprob(2,i,m) = 1.0; + end +else + obsprob(2,:,:) = mk_stochastic(rand(motif_length, Osize)); +end + +if 0 + Oargs = {'CPT', obsprob}; +else + % We use a minent prior for the emission distribution for the states in the motif model + % (but not the background model). This encourages nearly deterministic distributions. + % We create an index matrix (where M = motif length) + % [2 1 + % 2 2 + % ... + % 2 M] + % and then convert this to a list of integers, which + % specifies when to use the minent prior (Q1=2 specifies motif model). + M = motif_length; + ndx = [2*ones(M,1) (1:M)']; + pcases = subv2ind([2 motif_length], ndx); + Oargs = {'CPT', obsprob, 'prior_type', 'entropic', 'entropic_pcases', pcases}; +end + + + +[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ... + 'Oargs', Oargs, 'Ops', Qnodes(1:2), ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m new file mode 100644 index 00000000..b5822e10 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m @@ -0,0 +1,10 @@ +%bnet = mk_motif_hhmm('motif_pattern', 'acca', 'background', 't'); +bnet = mk_motif_hhmm('motif_pattern', 'accaggggga', 'background', []); + +chars = ['a', 'c', 'g', 't']; +Tmax = 100; + +for seqi=1:5 + evidence = cell2num(sample_dbn(bnet, 'length', Tmax)); + chars(evidence(end,:)) +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Entries new file mode 100644 index 00000000..6caae4a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Entries @@ -0,0 +1,8 @@ +/mk_abcd_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_arrow_alpha_hhmm3.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_hhmm2.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_hhmm3.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_hhmm3_args.m/1.1.1.1/Wed May 29 15:59:54 2002// +/motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/remove_hhmm_end_state.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Repository new file mode 100644 index 00000000..cc9acc63 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Old diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_abcd_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_abcd_hhmm.m new file mode 100644 index 00000000..330bc304 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_abcd_hhmm.m @@ -0,0 +1,109 @@ +% Make the HHMM in Figure 1 of the NIPS'01 paper + +Qsize = [2 3 2]; +D = 3; + +% transprob{d}(i,k,j), transprob{1}(i,j) +% termprob{d}(k,j), termprob{1}(1,j) +% startprob{d}(k,j), startprob{1}(1,j) +% obsprob(k, o) for discrete outputs + +% LEVEL 1 +% 1 2 e +A{1} = [0 0 1; + 0 0 1]; +[transprob{1}, termprob{1}] = remove_hhmm_end_state(A{1}); +startprob{1} = [0.5 0.5]; +Q1args = {'startprob', startprob{1}, 'transprob', transprob{1}}; + +% LEVEL 2 +A{2} = zeros(Qsize(2), Qsize(1), Qsize(2)+1); + +% 1 2 3 e +A{2}(:,1,:) = [0 1 0 0 + 0 0 1 0 + 0 0 0 1]; + +% 1 2 3 e +A{2}(:,2,:) = [0 1 0 0 + 0 0 1 0 + 0 0 0 1]; + +[transprob{2}, termprob{2}] = remove_hhmm_end_state(A{2}); + +% always enter level 2 in state 1 +startprob{2} = [1 0 0 + 1 0 0]; + +Q2args = {'startprob', startprob{2}, 'transprob', transprob{2}}; +F2args = {'CPT', termprob{2}}; + + +% LEVEL 3 + +A{3} = zeros([Qsize(3) Qsize(1:2) Qsize(3)+1]); +endstate = Qsize(3)+1; +% Qt-1(3) Qt(1) Qt(2) Qt(3) +% 1 2 e +A{3}(1, 1, 1, endstate) = 1.0; +A{3}(:, 1, 2, :) = [0.0 1.0 0.0 + 0.5 0.0 0.5]; +A{3}(1, 1, 3, endstate) = 1.0; + +A{3}(1, 2, 1, endstate) = 1.0; +A{3}(:, 2, 2, :) = [0.0 1.0 0.0 + 0.5 0.0 0.5]; +A{3}(1, 2, 3, endstate) = 1.0; + +A{3} = reshape(A{3}, [Qsize(3) prod(Qsize(1:2)) Qsize(3)+1]); +[transprob{3}, termprob{3}] = remove_hhmm_end_state(A{3}); + +% define the vertical entry points to level 3 +startprob{3} = zeros(Qsize); +% Q1 Q2 Q3 +startprob{3}(1, 1, 1) = 1.0; +startprob{3}(1, 2, 1) = 1.0; +startprob{3}(1, 3, 1) = 1.0; + +startprob{3}(2, 1, 1) = 1.0; +startprob{3}(2, 2, 1) = 1.0; +startprob{3}(2, 3, 1) = 1.0; + +startprob{3} = reshape(startprob{3}, prod(Qsize(1:2)), Qsize(3)); + +chars = ['a', 'b', 'c', 'd', 'x', 'y']; +Osize = length(chars); + +obsprob = zeros([Qsize Osize]); +% 1 2 3 O +obsprob(1,1,1,find(chars == 'a')) = 1.0; + +obsprob(1,2,1,find(chars == 'x')) = 1.0; +obsprob(1,2,2,find(chars == 'y')) = 1.0; + +obsprob(1,3,1,find(chars == 'b')) = 1.0; + +obsprob(2,1,1,find(chars == 'c')) = 1.0; + +obsprob(2,2,1,find(chars == 'x')) = 1.0; +obsprob(2,2,2,find(chars == 'y')) = 1.0; + +obsprob(2,3,1,find(chars == 'd')) = 1.0; + +obsprob = reshape(obsprob, prod(Qsize), Osize); + +[intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D); + +hhmm.Qnodes = Qnodes; +hhmm.Fnodes = Fnodes; +hhmm.Onode = Onode; +hhmm.D = D; +hhmm.Qsize = Qsize; +hhmm.Osize = Osize; +hhmm.startprob = startprob; +hhmm.transprob = transprob; +hhmm.termprob = termprob; +hhmm.obsprob = obsprob; +hhmm.A = A; + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_arrow_alpha_hhmm3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_arrow_alpha_hhmm3.m new file mode 100644 index 00000000..ba1aa8cb --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_arrow_alpha_hhmm3.m @@ -0,0 +1,86 @@ +% Make the following HHMM +% +% LH RH +% / \ +% / \ +% LR -> UD -> RL -> DU RL -> UD -> LR -> DU +% \ +% \ +% Q1 -> Q2 +% +% where level 1 is fully interconnected (not shown) +% level 2 is left-right +% and each model at level 3 is a 2 state LR shared HMM + +Qsizes = [2 4 2]; +D = 3; + +% LEVEL 1 + +startprob1 = 'ergodic'; +transprob1 = 'ergodic'; + + +% LEVEL 2 + +startprob = zeros(2, 4); +% Q1 Q2 +startprob(1, 1) = 1; +startprob(2, 3) = 1; + +transprob = zeros(2, 4, 4); +transprob(1,:,:) = [0 1 0 0 + 0 0 1 0 + 0 0 0 1 + 0 0 0 1]; +transprob(2,:,:) = [0 0 0 1 + 1 0 0 0 + 0 1 0 0 + 0 0 0 1]; + +Q2args = {'startprob', startprob, 'transprob', transprob}; + +% always terminate in state 4 (default) +% F2args + +% LEVEL 3 + +% Defaults are fine: always start in state 1, left-right model, finish in state 2 + + +% OBS LEVEl + +chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; +Osize = length(chars); + +obsprob = zeros([4 2 Osize]); +% Q2 Q3 O +obsprob(1, 1, find(chars == 'L')) = 1.0; +obsprob(1, 2, find(chars == 'l')) = 1.0; + +obsprob(2, 1, find(chars == 'U')) = 1.0; +obsprob(2, 2, find(chars == 'u')) = 1.0; + +obsprob(3, 1, find(chars == 'R')) = 1.0; +obsprob(3, 2, find(chars == 'r')) = 1.0; + +obsprob(4, 1, find(chars == 'D')) = 1.0; +obsprob(4, 2, find(chars == 'd')) = 1.0; + +Oargs = {'CPT', obsprob}; + + +bnet = mk_hhmm3('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', 1, 'Oargs', Oargs, 'Q1args', Q1args, 'Q2args', Q2args); + +T = 20; +usecell = 0; +evidence = sample_dbn(bnet, T, usecell); +%chars(evidence(end,:)) + +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; + +pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, obs, chars); + +eclass = bnet.equiv_class; +S=struct(bnet.CPD{eclass(Q2,2)}) diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm2.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm2.m new file mode 100644 index 00000000..032015ba --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm2.m @@ -0,0 +1,111 @@ +function bnet = mk_hhmm2(varargin) +% MK_HHMM2 Make a 2 level Hierarchical HMM +% bnet = mk_hhmm2(...) +% +% 2-layer hierarchical HMM (node numbers in parens) +% +% Q1(1) ---------> Q1(5) +% / | \ / | +% | | v / | +% | | F2(3) --- / | +% | | ^ \ | +% | | / \ | +% | v \ v +% | Q2(2)--------> Q2 (6) +% | | +% \ | +% v v +% O(4) +% +% +% Optional arguments [default] +% +% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0] +% obsCPT - CPT(o,q1,q2) params for O ['rnd'] +% mu - mu(:,q1,q2) params for O [ [] ] +% Sigma - Sigma(:,q1,q2) params for O [ [] ] +% +% F2toQ1 - 1 if Q2 is an hhmm_CPD, 0 if F2 -> Q2 arc is absent, so level 2 never resets [1] +% Q1args - arguments to be passed to the constructors for Q1(t=2) [ {} ] +% Q2args - arguments to be passed to the constructors for Q2(t=2) [ {} ] +% +% F2 only turns on (wp 0.5) when Q2 enters its final state. +% Q1 (slice 1) is clamped to be uniform. +% Q2 (slice 1) is clamped to always start in state 1. + +[os nmodels nstates] = size(mu); + +ss = 4; +Q1 = 1; Q2 = 2; F2 = 3; obs = 4; +Qnodes = [Q1 Q2]; +names = {'Q1', 'Q2', 'F2', 'obs'}; +intra = zeros(ss); +intra(Q1, [Q2 F2 obs]) = 1; +intra(Q2, [F2 obs]) = 1; + +inter = zeros(ss); +inter(Q1,Q1) = 1; +inter(F2,Q1) = 1; +if F2toQ2 + inter(F2,Q2)=1; +end +inter(Q2,Q2) = 1; + +ns = zeros(1,ss); + +ns(Q1) = nmodels; +ns(Q2) = nstates; +ns(F2) = 2; +ns(obs) = os; + +dnodes = [Q1 Q2 F2]; +if discrete_obs + dnodes = [dnodes obs]; +end +onodes = [obs]; + +bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names); +eclass = bnet.equiv_class; + +% SLICE 1 + +% We clamp untied nodes in the first slice, since their params can't be estimated +% from just one sequence + +% uniform prior on initial model +CPT = normalise(ones(1,nmodels)); +bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', CPT, 'adjustable', 0); + +% each model always starts in state 1 +CPT = zeros(ns(Q1), ns(Q2)); +CPT(:, 1) = 1.0; +bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', CPT, 'adjustable', 0); + +% Termination probability +CPT = zeros(ns(Q1), ns(Q2), 2); +if 1 + % Each model can only terminate in its final state. + % 0 params will remain 0 during EM, thus enforcing this constraint. + CPT(:, :, 1) = 1.0; % all states turn F off ... + p = 0.5; + CPT(:, ns(Q2), 2) = p; % except the last one + CPT(:, ns(Q2), 1) = 1-p; +end +bnet.CPD{eclass(F2,1)} = tabular_CPD(bnet, F2, 'CPT', CPT); + +if discrete_obs + bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, obs_args{:}); +else + bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, obs_args{:}); +end + +% SLICE 2 + + +bnet.CPD{eclass(Q1,2)} = hhmm_CPD(bnet, Q1+ss, Qnodes, 1, D, 'args', Q1args); + +if F2toQ2 + bnet.CPD{eclass(Q2,2)} = hhmmQD_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:}); +else + bnet.CPD{eclass(Q2,2)} = tabular_CPD(bnet, Q2+ss, Q2args{:}); +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3.m new file mode 100644 index 00000000..5b2cd6aa --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3.m @@ -0,0 +1,181 @@ +function bnet = mk_hhmm3(varargin) +% MK_HHMM3 Make a 3 level Hierarchical HMM +% bnet = mk_hhmm3(...) +% +% 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs. +% This enforces sub-models (which differ only in their Q1 index) to be shared. +% Also, we enforce the fact that each model always starts in its initial state +% and only finishes in its final state. However, the prob. of finishing (as opposed to +% self-transitioning to the final state) can be learned. +% The fact that we always finish from the same state means we do not need to condition +% F(i) on Q(i-1), since finishing prob is indep of calling context. +% +% The DBN is the same as Fig 10 in my tech report. +% +% Q1 ----------> Q1 +% | / | +% | / | +% | F2 ------- | +% | ^ \ | +% | /| \ | +% v | v v +% Q2-| --------> Q2 +% /| | ^ +% / | | /| +% | | F3 ---------/ | +% | | ^ \ | +% | v / v +% | Q3 -----------> Q3 +% | | +% \ | +% v v +% O +% +% +% Optional arguments in name/value format [default] +% +% Qsizes - sizes at each level [ none ] +% Osize - size of O node [ none ] +% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0] +% Oargs - cell array of args to pass to the O CPD [ {} ] +% transprob1 - transprob1(i,j) = P(Q1(t)=j|Q1(t-1)=i) ['ergodic'] +% startprob1 - startprob1(j) = P(Q1(t)=j) ['leftstart'] +% transprob2 - transprob2(i,k,j) = P(Q2(t)=j|Q2(t-1)=i,Q1(t)=k) ['leftright'] +% startprob2 - startprob2(k,j) = P(Q2(t)=j|Q1(t)=k) ['leftstart'] +% termprob2 - termprob2(j,f) = P(F2(t)=f|Q2(t)=j) ['rightstop'] +% transprob3 - transprob3(i,k,j) = P(Q3(t)=j|Q3(t-1)=i,Q2(t)=k) ['leftright'] +% startprob3 - startprob3(k,j) = P(Q3(t)=j|Q2(t)=k) ['leftstart'] +% termprob3 - termprob3(j,f) = P(F3(t)=f|Q3(t)=j) ['rightstop'] +% +% leftstart means the model always starts in state 1. +% rightstop means the model always finished in its last state (Qsize(d)). +% +% Q1:Q3 in slice 1 are of type tabular_CPD +% Q1:Q3 in slice 2 are of type hhmmQ_CPD. +% F2 is of type hhmmF_CPD, F3 is of type tabular_CPD. + +ss = 6; D = 3; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; +names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'}; + +intra = zeros(ss); +intra(Q1, Q2) = 1; +intra(Q2, [F2 Q3 obs]) = 1; +intra(Q3, [F3 obs]) = 1; +intra(F3, F2) = 1; + +inter = zeros(ss); +inter(Q1,Q1) = 1; +inter(Q2,Q2) = 1; +inter(Q3,Q3) = 1; +inter(F2,[Q1 Q2]) = 1; +inter(F3,[Q2 Q3]) = 1; + + +% get sizes of nodes +args = varargin; +nargs = length(args); +Qsizes = []; +Osize = 0; +for i=1:2:nargs + switch args{i}, + case 'Qsizes', Qsizes = args{i+1}; + case 'Osize', Osize = args{i+1}; + end +end +if isempty(Qsizes), error('must specify Qsizes'); end +if Osize==0, error('must specify Osize'); end + +% set default params +discrete_obs = 0; +Oargs = {}; +startprob1 = 'ergodic'; +startprob2 = 'leftstart'; +startprob3 = 'leftstart'; +transprob1 = 'ergodic'; +transprob2 = 'leftright'; +transprob3 = 'leftright'; +termprob2 = 'rightstop'; +termprob3 = 'rightstop'; + + +for i=1:2:nargs + switch args{i}, + case 'discrete_obs', discrete_obs = args{i+1}; + case 'Oargs', Oargs = args{i+1}; + case 'Q1args', Q1args = args{i+1}; + case 'Q2args', Q2args = args{i+1}; + case 'Q3args', Q3args = args{i+1}; + case 'F2args', F2args = args{i+1}; + case 'F3args', F3args = args{i+1}; + end +end + + +ns = zeros(1,ss); +ns(Qnodes) = Qsizes; +ns(obs) = Osize; +ns(Fnodes) = 2; + +dnodes = [Qnodes Fnodes]; +if discrete_obs + dnodes = [dnodes obs]; +end +onodes = [obs]; + +bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names); +eclass = bnet.equiv_class; + +if strcmp(startprob1, 'ergodic') + startprob1 = normalise(ones(1,ns(Q1))); +end +if strcmp(startprob2, 'leftstart') + startprob2 = zeros(ns(Q1), ns(Q2)); + starpbrob2(:, 1) = 1.0; +end +if strcmp(startprob3, 'leftstart') + startprob3 = zeros(ns(Q2), ns(Q3)); + starpbrob3(:, 1) = 1.0; +end + +if strcmp(termprob2, 'rightstop') + p = 0.9; + termprob2 = zeros(Qsize(2),2); + termprob2(:, 2) = p; + termprob2(:, 1) = 1-p; + termprob2(1:(Qsize(2)-1), 1) = 1; +end +if strcmp(termprob3, 'rightstop') + p = 0.9; + termprob3 = zeros(Qsize(3),2); + termprob3(:, 2) = p; + termprob3(:, 1) = 1-p; + termprob3(1:(Qsize(3)-1), 1) = 1; +end + + +% SLICE 1 + +% We clamp untied nodes in the first slice, since their params can't be estimated +% from just one sequence + +bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', startprob1, 'adjustable', 0); +bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', startprob2, 'adjustable', 0); +bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', startprob3, 'adjustable', 0); + +bnet.CPD{eclass(F2,1)} = hhmmF_CPD(bnet, F2, Qnodes, 2, D, 'termprob', termprob2); +bnet.CPD{eclass(F3,1)} = tabular_CPD(bnet, F3, 'CPT', termprob3); + +if discrete_obs + bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:}); +else + bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:}); +end + +% SLICE 2 + +bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, 'transprob', transprob1, 'startprob', startprob1); +bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, 'transprob', transprob2, 'startprob', startprob2); +bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, 'transprob', transprob3, 'startprob', startprob3); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3_args.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3_args.m new file mode 100644 index 00000000..bc3ec886 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/mk_hhmm3_args.m @@ -0,0 +1,165 @@ +function bnet = mk_hhmm3(varargin) +% MK_HHMM3 Make a 3 level Hierarchical HMM +% bnet = mk_hhmm3(...) +% +% 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs. +% This enforces sub-models (which differ only in their Q1 index) to be shared. +% Also, we enforce the fact that each model always starts in its initial state +% and only finishes in its final state. However, the prob. of finishing (as opposed to +% self-transitioning to the final state) can be learned. +% The fact that we always finish from the same state means we do not need to condition +% F(i) on Q(i-1), since finishing prob is indep of calling context. +% +% The DBN is the same as Fig 10 in my tech report. +% +% Q1 ----------> Q1 +% | / | +% | / | +% | F2 ------- | +% | ^ \ | +% | /| \ | +% v | v v +% Q2-| --------> Q2 +% /| | ^ +% / | | /| +% | | F3 ---------/ | +% | | ^ \ | +% | v / v +% | Q3 -----------> Q3 +% | | +% \ | +% v v +% O +% +% Q1 (slice 1) is clamped to be uniform. +% Q2 (slice 1) is clamped to always start in state 1. +% Q3 (slice 1) is clamped to always start in state 1. +% F3 by default will only finish if Q3 is in its last state (F3 is a tabular_CPD) +% F2 by default gets the default hhmmF_CPD params. +% Q1:Q3 (slice 2) by default gets the default hhmmQ_CPD params. +% O by default gets the default tabular/Gaussian params. +% +% Optional arguments in name/value format [default] +% +% Qsizes - sizes at each level [ none ] +% Osize - size of O node [ none ] +% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0] +% Oargs - cell array of args to pass to the O CPD [ {} ] +% Q1args - args to be passed to constructor for Q1 (slice 2) [ {} ] +% Q2args - args to be passed to constructor for Q2 (slice 2) [ {} ] +% Q3args - args to be passed to constructor for Q3 (slice 2) [ {} ] +% F2args - args to be passed to constructor for F2 [ {} ] +% F3args - args to be passed to constructor for F3 [ {'CPT', finish in last Q3 state} ] +% + +ss = 6; D = 3; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; +names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'}; + +intra = zeros(ss); +intra(Q1, Q2) = 1; +intra(Q2, [F2 Q3 obs]) = 1; +intra(Q3, [F3 obs]) = 1; +intra(F3, F2) = 1; + +inter = zeros(ss); +inter(Q1,Q1) = 1; +inter(Q2,Q2) = 1; +inter(Q3,Q3) = 1; +inter(F2,[Q1 Q2]) = 1; +inter(F3,[Q2 Q3]) = 1; + + +% get sizes of nodes +args = varargin; +nargs = length(args); +Qsizes = []; +Osize = 0; +for i=1:2:nargs + switch args{i}, + case 'Qsizes', Qsizes = args{i+1}; + case 'Osize', Osize = args{i+1}; + end +end +if isempty(Qsizes), error('must specify Qsizes'); end +if Osize==0, error('must specify Osize'); end + +% set default params +discrete_obs = 0; +Oargs = {}; +Q1args = {}; +Q2args = {}; +Q3args = {}; +F2args = {}; + +% P(Q3, F3) +CPT = zeros(Qsizes(3), 2); +% Each model can only terminate in its final state. +% 0 params will remain 0 during EM, thus enforcing this constraint. +CPT(:, 1) = 1.0; % all states turn F off ... +p = 0.5; +CPT(Qsizes(3), 2) = p; % except the last one +CPT(Qsizes(3), 1) = 1-p; +F3args = {'CPT', CPT}; + +for i=1:2:nargs + switch args{i}, + case 'discrete_obs', discrete_obs = args{i+1}; + case 'Oargs', Oargs = args{i+1}; + case 'Q1args', Q1args = args{i+1}; + case 'Q2args', Q2args = args{i+1}; + case 'Q3args', Q3args = args{i+1}; + case 'F2args', F2args = args{i+1}; + case 'F3args', F3args = args{i+1}; + end +end + +ns = zeros(1,ss); +ns(Qnodes) = Qsizes; +ns(obs) = Osize; +ns(Fnodes) = 2; + +dnodes = [Qnodes Fnodes]; +if discrete_obs + dnodes = [dnodes obs]; +end +onodes = [obs]; + +bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names); +eclass = bnet.equiv_class; + +% SLICE 1 + +% We clamp untied nodes in the first slice, since their params can't be estimated +% from just one sequence + +% uniform prior on initial model +CPT = normalise(ones(1,ns(Q1))); +bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', CPT, 'adjustable', 0); + +% each model always starts in state 1 +CPT = zeros(ns(Q1), ns(Q2)); +CPT(:, 1) = 1.0; +bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', CPT, 'adjustable', 0); + +% each model always starts in state 1 +CPT = zeros(ns(Q2), ns(Q3)); +CPT(:, 1) = 1.0; +bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', CPT, 'adjustable', 0); + +bnet.CPD{eclass(F2,1)} = hhmmF_CPD(bnet, F2, Qnodes, 2, D, F2args{:}); + +bnet.CPD{eclass(F3,1)} = tabular_CPD(bnet, F3, F3args{:}); + +if discrete_obs + bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:}); +else + bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:}); +end + +% SLICE 2 + +bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, Q1args{:}); +bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:}); +bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, Q3args{:}); diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/motif_hhmm.m new file mode 100644 index 00000000..10144b58 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/motif_hhmm.m @@ -0,0 +1,95 @@ +% Make the following HHMM +% +% S1 <----------------------> S2 +% | | +% | | +% M1 -> M2 -> M3 -> end B1 -> end +% +% where Mi represents the i'th letter in the motif +% and B is the background state. +% Si chooses between running the motif or the background. +% The Si and B states have self loops (not shown). + +if 0 +seed = 0; +rand('state', seed); +randn('state', seed); +end + +chars = ['a', 'c', 'g', 't']; +Osize = length(chars); + +motif_pattern = 'acca'; +motif_length = length(motif_pattern); +Qsize = [2 motif_length]; +Qnodes = 1:2; +D = 2; +transprob = cell(1,D); +termprob = cell(1,D); +startprob = cell(1,D); + +% startprob{d}(k,j), startprob{1}(1,j) +% transprob{d}(i,k,j), transprob{1}(i,j) +% termprob{d}(k,j) + + +% LEVEL 1 + +startprob{1} = zeros(1, 2); +startprob{1} = [1 0]; % always start in the background model + +% When in the background state, we stay there with high prob +% When in the motif state, we immediately return to the background state. +transprob{1} = [0.8 0.2; + 1.0 0.0]; + + +% LEVEL 2 +startprob{2} = 'leftstart'; % both submodels start in substate 1 +transprob{2} = zeros(motif_length, 2, motif_length); +termprob{2} = zeros(2, motif_length); + +% In the background model, we only use state 1. +transprob{2}(1,1,1) = 1; % self loop +termprob{2}(1,1) = 0.2; % prob transition to end state + +% Motif model +transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); +termprob{2}(2,end) = 1.0; % last state immediately terminates + + +% OBS LEVEl + +obsprob = zeros([Qsize Osize]); +if 0 + % uniform background model + obsprob(1,1,:) = normalise(ones(Osize,1)); +else + % deterministic background model (easy to see!) + m = find(chars=='t'); + obsprob(1,1,m) = 1.0; +end +if 1 + % initialise with true motif (cheating) + for i=1:motif_length + m = find(chars == motif_pattern(i)); + obsprob(2,i,m) = 1.0; + end +end + +Oargs = {'CPT', obsprob}; + +[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ... + 'Oargs', Oargs, 'Ops', Qnodes(1:2), ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + + +Tmax = 20; +usecell = 0; + +for seqi=1:5 + evidence = sample_dbn(bnet, Tmax, usecell); + chars(evidence(end,:)) + %T = size(evidence, 2) + %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars); +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/remove_hhmm_end_state.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/remove_hhmm_end_state.m new file mode 100644 index 00000000..2bf10d72 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Old/remove_hhmm_end_state.m @@ -0,0 +1,37 @@ +function [transprob, termprob] = remove_hhmm_end_state(A) +% REMOVE_END_STATE Infer transition and termination probabilities from automaton with an end state +% [transprob, termprob] = remove_end_state(A) +% A(i,k,j) = Pr( i->j | Qps=k), where i in 1:Q, j in 1:(Q+1), and Q+1 is the end state + +if ndims(A)==2 % top level + Q = size(A,1); + transprob = A(:,1:Q); + termprob = A(:,Q+1)'; + + % rescale + for i=1:Q + for j=1:Q + denom = (1-termprob(i)); + denom = denom + (denom==0)*eps; + transprob(i,j) = transprob(i,j) / denom; + end + end +else + Q = size(A,1); + Qk = size(A,2); + transprob = A(:, :, 1:Q); + termprob = A(:,:,Q+1)'; + + % rescale + for k=1:Qk + for i=1:Q + for j=1:Q + denom = (1-termprob(k,i)); + denom = denom + (denom==0)*eps; + transprob(i,k,j) = transprob(i,k,j) / denom; + end + end + end + +end + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Entries new file mode 100644 index 00000000..fdee19ae --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Entries @@ -0,0 +1,14 @@ +/get_square_data.m/1.1.1.1/Wed May 29 15:59:54 2002// +/hhmm_inference.m/1.1.1.1/Wed May 29 15:59:54 2002// +/is_F2_true_D3.m/1.1.1.1/Wed May 29 15:59:54 2002// +/learn_square_hhmm_cts.m/1.1.1.1/Thu Jun 20 00:19:22 2002// +/learn_square_hhmm_discrete.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/plot_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_square_hhmm_cts.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_square_hhmm_discrete.m/1.1.1.1/Wed May 29 15:59:54 2002// +/square4.mat/1.1.1.1/Wed May 29 15:59:54 2002// +/square4_cases.mat/1.1.1.1/Wed May 29 15:59:54 2002// +/test_square_fig.m/1.1.1.1/Wed May 29 15:59:54 2002// +/test_square_fig.mat/1.1.1.1/Wed May 29 15:59:54 2002// +D/Old//// diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Repository new file mode 100644 index 00000000..e926a0d5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Square diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Entries new file mode 100644 index 00000000..6d415d0d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Entries @@ -0,0 +1,5 @@ +/learn_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/mk_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/plot_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_square_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Repository new file mode 100644 index 00000000..47df1a8c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Square/Old diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/learn_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/learn_square_hhmm.m new file mode 100644 index 00000000..695ae047 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/learn_square_hhmm.m @@ -0,0 +1,294 @@ +% Learn a 3 level HHMM similar to mk_square_hhmm + +% Because startprob should be shared for t=1:T, +% but in the DBN is shared for t=2:T, we train using a single long sequence. + +discrete_obs = 0; +supervised = 1; +obs_finalF2 = 0; +% It is not possible to observe F2 if we learn +% because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume +% the F nodes are always hidden (for speed). +% However, for generating, we might want to set the final F2=true +% to force all subroutines to finish. + +ss = 6; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; + +seed = 1; +rand('state', seed); +randn('state', seed); + +if discrete_obs + Qsizes = [2 4 2]; +else + Qsizes = [2 4 1]; +end + +D = 3; +Qnodes = 1:D; +startprob = cell(1,D); +transprob = cell(1,D); +termprob = cell(1,D); + +startprob{1} = 'unif'; +transprob{1} = 'unif'; + +% In the unsupervised case, it is essential that we break symmetry +% in the initial param estimates. +%startprob{2} = 'unif'; +%transprob{2} = 'unif'; +%termprob{2} = 'unif'; +startprob{2} = 'rnd'; +transprob{2} = 'rnd'; +termprob{2} = 'rnd'; + +leftright = 0; +if leftright + % Initialise base-level models as left-right. + % If we initialise with delta functions, + % they will remain delat funcitons after learning + startprob{3} = 'leftstart'; + transprob{3} = 'leftright'; + termprob{3} = 'rightstop'; +else + % If we want to be able to run a base-level model backwards... + startprob{3} = 'rnd'; + transprob{3} = 'rnd'; + termprob{3} = 'rnd'; +end + +if discrete_obs + % Initialise observations of lowest level primitives in a way which we can interpret + chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; + L=find(chars=='L'); l=find(chars=='l'); + U=find(chars=='U'); u=find(chars=='u'); + R=find(chars=='R'); r=find(chars=='r'); + D=find(chars=='D'); d=find(chars=='d'); + Osize = length(chars); + + p = 0.9; + obsprob = (1-p)*ones([4 2 Osize]); + % Q2 Q3 O + obsprob(1, 1, L) = p; + obsprob(1, 2, l) = p; + obsprob(2, 1, U) = p; + obsprob(2, 2, u) = p; + obsprob(3, 1, R) = p; + obsprob(3, 2, r) = p; + obsprob(4, 1, D) = p; + obsprob(4, 2, d) = p; + obsprob = mk_stochastic(obsprob); + Oargs = {'CPT', obsprob}; + +else + % Initialise means of lowest level primitives in a way which we can interpret + % These means are little vectors in the east, south, west, north directions. + % (left-right=east, up-down=south, right-left=west, down-up=north) + Osize = 2; + mu = zeros(2, Qsizes(2), Qsizes(3)); + noise = 0; + scale = 3; + for q3=1:Qsizes(3) + mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1); + end + Sigma = repmat(reshape(scale*eye(2), [2 2 1 1 ]), [1 1 Qsizes(2) Qsizes(3)]); + Oargs = {'mean', mu, 'cov', Sigma, 'cov_type', 'diag'}; +end + +bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs,... + 'Oargs', Oargs, 'Ops', Qnodes(2:3), ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + +if supervised + bnet.observed = [Q1 Q2 Onode]; +else + bnet.observed = [Onode]; +end + +if obs_finalF2 + engine = jtree_dbn_inf_engine(bnet); + % can't use ndx version because sometimes F2 is hidden, sometimes observed + error('can''t observe F when learning') +else + if supervised + engine = jtree_ndx_dbn_inf_engine(bnet); + else + engine = jtree_hmm_inf_engine(bnet); + end +end + +if discrete_obs + % generate some synthetic data (easier to debug) + cases = {}; + + T = 8; + ev = cell(ss, T); + ev(Onode,:) = num2cell([L l U u R r D d]); + if supervised + ev(Q1,:) = num2cell(1*ones(1,T)); + ev(Q2,:) = num2cell( [1 1 2 2 3 3 4 4]); + end + cases{1} = ev; + cases{3} = ev; + + T = 8; + ev = cell(ss, T); + if leftright % base model is left-right + ev(Onode,:) = num2cell([R r U u L l D d]); + else + ev(Onode,:) = num2cell([r R u U l L d D]); + end + if supervised + ev(Q1,:) = num2cell(2*ones(1,T)); + ev(Q2,:) = num2cell( [3 3 2 2 1 1 4 4]); + end + + cases{2} = ev; + cases{4} = ev; + + if obs_finalF2 + for i=1:length(cases) + T = size(cases{i},2); + cases{i}(F2,T)={2}; % force F2 to be finished at end of seq + end + end + + if 0 + ev = cases{4}; + engine2 = enter_evidence(engine2, ev); + T = size(ev,2); + for t=1:T + m=marginal_family(engine2, F2, t); + fprintf('t=%d\n', t); + reshape(m.T, [2 2]) + end + end + + % [bnet2, LL] = learn_params_dbn_em(engine, cases, 'max_iter', 10); + long_seq = cat(2, cases{:}); + [bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 200); + + % figure out which subsequence each model is responsible for + mpe = calc_mpe_dbn(engine2, long_seq); + pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, chars); + +else + load 'square4_cases' % cases{seq}{i,t} for i=1:ss + %plot_square_hhmm(cases{1}) + %long_seq = cat(2, cases{:}); + train_cases = cases(1:2); + long_seq = cat(2, train_cases{:}); + if ~supervised + T = size(long_seq,2); + for t=1:T + long_seq{Q1,t} = []; + long_seq{Q2,t} = []; + end + end + [bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 100); + + CPDO=struct(bnet2.CPD{eclass(Onode,1)}); + mu = CPDO.mean; + Sigma = CPDO.cov; + CPDO_full = CPDO; + + % force diagonal covs after training + for k=1:size(Sigma,3) + Sigma(:,:,k) = diag(diag(Sigma(:,:,k))); + end + bnet2.CPD{6} = set_fields(bnet.CPD{6}, 'cov', Sigma); + + if 0 + % visualize each model by concatenating means for each model for nsteps in a row + nsteps = 5; + ev = cell(ss, nsteps*prod(Qsizes(2:3))); + t = 1; + for q2=1:Qsizes(2) + for q3=1:Qsizes(3) + for i=1:nsteps + ev{Onode,t} = mu(:,q2,q3); + ev{Q2,t} = q2; + t = t + 1; + end + end + end + plot_square_hhmm(ev) + end + + % bnet3 is the same as the learned model, except we will use it in testing mode + if supervised + bnet3 = bnet2; + bnet3.observed = [Onode]; + engine3 = hmm_inf_engine(bnet3); + %engine3 = jtree_ndx_dbn_inf_engine(bnet3); + else + bnet3 = bnet2; + engine3 = engine2; + end + + if 0 + % segment whole sequence + mpe = calc_mpe_dbn(engine3, long_seq); + pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []); + end + + % segment each sequence + test_cases = cases(3:4); + for i=1:2 + ev = test_cases{i}; + T = size(ev, 2); + for t=1:T + ev{Q1,t} = []; + ev{Q2,t} = []; + end + mpe = calc_mpe_dbn(engine3, ev); + subplot(1,2,i) + plot_square_hhmm(mpe) + %pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []); + q1s = cell2num(mpe(Q1,:)); + h = hist(q1s, 1:Qsizes(1)); + map_q1 = argmax(h); + str = sprintf('test seq %d is of type %d\n', i, map_q1); + title(str) + end + +end + +if 0 +% Estimate gotten by couting transitions in the labelled data +% Note that a self transition shouldnt count if F2=off. +Q2ev = cell2num(ev(Q2,:)); +Q2a = Q2ev(1:end-1); +Q2b = Q2ev(2:end); +counts = compute_counts([Q2a; Q2b], [4 4]); +end + +eclass = bnet2.equiv_class; +CPDQ1=struct(bnet2.CPD{eclass(Q1,2)}); +CPDQ2=struct(bnet2.CPD{eclass(Q2,2)}); +CPDQ3=struct(bnet2.CPD{eclass(Q3,2)}); +CPDF2=struct(bnet2.CPD{eclass(F2,1)}); +CPDF3=struct(bnet2.CPD{eclass(F3,1)}); + + +A=add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2)); +squeeze(A(:,1,:)) +squeeze(A(:,2,:)) +CPDQ2.startprob + +if 0 +S=struct(CPDF2.sub_CPD_term); +S.nsamples +reshape(S.counts, [2 4 2]) +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/mk_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/mk_square_hhmm.m new file mode 100644 index 00000000..608b6784 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/mk_square_hhmm.m @@ -0,0 +1,183 @@ +function bnet = mk_square_hhmm(discrete_obs, true_params, topright) + +% Make a 3 level HHMM described by the following grammar +% +% Square -> CLK | CCK % clockwise or counterclockwise +% CLK -> LR UD RL DU start on top left (1 2 3 4) +% CCK -> RL UD LR DU if start at top right (3 2 1 4) +% CCK -> UD LR DU RL if start at top left (2 1 4 3) +% +% LR = left-right, UD = up-down, RL = right-left, DU = down-up +% LR, UD, RL, DU are sub HMMs. +% +% For discrete observations, the subHMMs are 2-state left-right. +% LR emits L then l, etc. +% +% For cts observations, the subHMMs are 1 state. +% LR emits a vector in the -> direction, with a little noise. +% Since there is no constraint that we remain in the LR state as long as the RL state, +% the sides of the square might have different lengths, +% so the result is not really a square! +% +% If true_params = 0, we use random parameters at the top 2 levels +% (ready for learning). At the bottom level, we use noisy versions +% of the "true" observations. +% +% If topright=1, counter-clockwise starts at top right, not top left +% This example was inspired by Ivanov and Bobick. + +if nargin < 3, topright = 1; end + +if 1 % discrete_obs + Qsizes = [2 4 2]; +else + Qsizes = [2 4 1]; +end + +D = 3; +Qnodes = 1:D; +startprob = cell(1,D); +transprob = cell(1,D); +termprob = cell(1,D); + +% LEVEL 1 + +startprob{1} = 'unif'; +transprob{1} = 'unif'; + +% LEVEL 2 + +if true_params + startprob{2} = zeros(2, 4); + startprob{2}(1, :) = [1 0 0 0]; + if topright + startprob{2}(2, :) = [0 0 1 0]; + else + startprob{2}(2, :) = [0 1 0 0]; + end + + transprob{2} = zeros(4, 2, 4); + + transprob{2}(:,1,:) = [0 1 0 0 + 0 0 1 0 + 0 0 0 1 + 0 0 0 1]; % 4->e + if topright + transprob{2}(:,2,:) = [0 0 0 1 + 1 0 0 0 + 0 1 0 0 + 0 0 0 1]; % 4->e + else + transprob{2}(:,2,:) = [0 0 0 1 + 1 0 0 0 + 0 0 1 0 % 3->e + 0 0 1 0]; + end + + %termprob{2} = 'rightstop'; + termprob{2} = zeros(2,4,2); + pfin = 0.8; + termprob{2}(1,:,2) = [0 0 0 pfin]; % finish in state 4 (DU) + termprob{2}(1,:,1) = 1 - [0 0 0 pfin]; + if topright + termprob{2}(2,:,2) = [0 0 0 pfin]; + termprob{2}(2,:,1) = 1 - [0 0 0 pfin]; + else + termprob{2}(2,:,2) = [0 0 pfin 0]; % finish in state 3 (RL) + termprob{2}(2,:,1) = 1 - [0 0 pfin 0]; + end +else + % In the unsupervised case, it is essential that we break symmetry + % in the initial param estimates. + %startprob{2} = 'unif'; + %transprob{2} = 'unif'; + %termprob{2} = 'unif'; + startprob{2} = 'rnd'; + transprob{2} = 'rnd'; + termprob{2} = 'rnd'; +end + +% LEVEL 3 + +if 1 | true_params + startprob{3} = 'leftstart'; + transprob{3} = 'leftright'; + termprob{3} = 'rightstop'; +else + % If we want to be able to run a base-level model backwards... + startprob{3} = 'rnd'; + transprob{3} = 'rnd'; + termprob{3} = 'rnd'; +end + + +% OBS LEVEl + +if discrete_obs + % Initialise observations of lowest level primitives in a way which we can interpret + chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; + L=find(chars=='L'); l=find(chars=='l'); + U=find(chars=='U'); u=find(chars=='u'); + R=find(chars=='R'); r=find(chars=='r'); + D=find(chars=='D'); d=find(chars=='d'); + Osize = length(chars); + + if true_params + p = 1; % makes each state fully observed + else + p = 0.9; + end + + obsprob = (1-p)*ones([4 2 Osize]); + % Q2 Q3 O + obsprob(1, 1, L) = p; + obsprob(1, 2, l) = p; + obsprob(2, 1, U) = p; + obsprob(2, 2, u) = p; + obsprob(3, 1, R) = p; + obsprob(3, 2, r) = p; + obsprob(4, 1, D) = p; + obsprob(4, 2, d) = p; + obsprob = mk_stochastic(obsprob); + Oargs = {'CPT', obsprob}; +else + % Initialise means of lowest level primitives in a way which we can interpret + % These means are little vectors in the east, south, west, north directions. + % (left-right=east, up-down=south, right-left=west, down-up=north) + Osize = 2; + mu = zeros(2, Qsizes(2), Qsizes(3)); + scale = 3; + if true_params + noise = 0; + else + noise = 0.5*scale; + end + for q3=1:Qsizes(3) + mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1); + end + Sigma = repmat(reshape(scale*eye(2), [2 2 1 1 ]), [1 1 Qsizes(2) Qsizes(3)]); + Oargs = {'mean', mu, 'cov', Sigma, 'cov_type', 'diag'}; +end + +if discrete_obs + selfprob = 0.5; +else + selfprob = 0.95; + % If less than this, it won't look like a square + % because it doesn't spend enough time in each state + % Unfortunately, the variance on durations (lengths of each side) + % is very large +end +bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ... + 'Oargs', Oargs, 'Ops', Qnodes(2:3), 'selfprob', selfprob, ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/plot_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/plot_square_hhmm.m new file mode 100644 index 00000000..e6701e45 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/plot_square_hhmm.m @@ -0,0 +1,32 @@ +function plot_square_hhmm(ev) +% Plot the square shape implicit in the evidence. +% ev{i,t} is the value of node i in slice t. +% The observed node contains a velocity (delta increment), which is converted +% into a position. +% The Q2 node specifies which model is used; each segment is color-coded +% in the order red, green, blue, black. + +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; + +delta = cell2num(ev(Onode,:)); % delta(:,t) +Q2label = cell2num(ev(Q2,:)); + +T = size(delta, 2); +pos = zeros(2,T+1); +clf +hold on +cols = {'r', 'g', 'b', 'k'}; +boundary = 0; +coli = 1; +for t=2:T+1 + pos(:,t) = pos(:,t-1) + delta(:,t-1); + plot(pos(1,t), pos(2,t), sprintf('%c.', cols{coli})); + if t < T + boundary = (Q2label(t) ~= Q2label(t-1)); + end + if boundary + coli = coli + 1; + coli = mod(coli-1, length(cols)) + 1; + end +end + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/sample_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/sample_square_hhmm.m new file mode 100644 index 00000000..a0f9007e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/Old/sample_square_hhmm.m @@ -0,0 +1,160 @@ + +seed = 0; +rand('state', seed); +randn('state', seed); + +discrete_obs = 1; +topright = 0; + +Qsizes = [2 4 2]; +D = 3; +Qnodes = 1:D; +startprob = cell(1,D); +transprob = cell(1,D); +termprob = cell(1,D); + +% LEVEL 1 + +startprob{1} = 'ergodic'; +transprob{1} = 'ergodic'; + +% LEVEL 2 + +startprob{2} = zeros(2, 4); +startprob{2}(1, :) = [1 0 0 0]; +if topright + startprob{2}(2, :) = [0 0 1 0]; +else + startprob{2}(2, :) = [0 1 0 0]; +end + +transprob{2} = zeros(4, 2, 4); + +transprob{2}(:,1,:) = [0 1 0 0 + 0 0 1 0 + 0 0 0 1 + 0 0 0 1]; % 4->e +if topright + transprob{2}(:,2,:) = [0 0 0 1 + 1 0 0 0 + 0 1 0 0 + 0 0 0 1]; % 4->e +else + transprob{2}(:,2,:) = [0 0 0 1 + 1 0 0 0 + 0 0 1 0 % 3->e + 0 0 1 0]; +end + +%termprob{2} = 'rightstop'; +termprob{2} = zeros(2,4,2); +pfin = 0.8; +termprob{2}(1,:,2) = [0 0 0 pfin]; % finish in state 4 (DU) +termprob{2}(1,:,1) = 1 - [0 0 0 pfin]; +if topright + termprob{2}(2,:,2) = [0 0 0 pfin]; + termprob{2}(2,:,1) = 1 - [0 0 0 pfin]; +else + termprob{2}(2,:,2) = [0 0 pfin 0]; % finish in state 3 (RL) + termprob{2}(2,:,1) = 1 - [0 0 pfin 0]; +end + +% LEVEL 3 + +startprob{3} = 'leftstart'; +transprob{3} = 'leftright'; +termprob{3} = 'rightstop'; + + +% OBS LEVEl + +if discrete_obs + chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; + L=find(chars=='L'); l=find(chars=='l'); + U=find(chars=='U'); u=find(chars=='u'); + R=find(chars=='R'); r=find(chars=='r'); + D=find(chars=='D'); d=find(chars=='d'); + Osize = length(chars); + + obsprob = zeros([4 2 Osize]); + % Q2 Q3 O + obsprob(1, 1, L) = 1.0; + obsprob(1, 2, l) = 1.0; + obsprob(2, 1, U) = 1.0; + obsprob(2, 2, u) = 1.0; + obsprob(3, 1, R) = 1.0; + obsprob(3, 2, r) = 1.0; + obsprob(4, 1, D) = 1.0; + obsprob(4, 2, d) = 1.0; + + Oargs = {'CPT', obsprob}; +else + Osize = 2; + mu = zeros(2, 4, 2); + noise = 0; + scale = 10; + for q3=1:2 + mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1); + end + for q3=1:2 + mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1); + end + for q3=1:2 + mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1); + end + for q3=1:2 + mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1); + end + Sigma = repmat(reshape(0.01*eye(2), [2 2 1 1 ]), [1 1 4 2]); + Oargs = {'mean', mu, 'cov', Sigma}; +end + +bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ... + 'Oargs', Oargs, 'Ops', Qnodes(2:3), ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + +if discrete_obs + Tmax = 30; +else + Tmax = 200; +end +usecell = ~discrete_obs; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; + +for seqi=1:3 + evidence = sample_dbn(bnet, Tmax, usecell, 'stop_sampling_F2'); + T = size(evidence, 2) + if discrete_obs + pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars); + else + pos = zeros(2,T+1); + delta = cell2num(evidence(Onode,:)); + clf + hold on + cols = {'r', 'g', 'k', 'b'}; + boundary = cell2num(evidence(F3,:))-1; + coli = 1; + for t=2:T+1 + pos(:,t) = pos(:,t-1) + delta(:,t-1); + plot(pos(1,t), pos(2,t), sprintf('%c.', cols{coli})); + if boundary(t-1) + coli = coli + 1; + coli = mod(coli-1, length(cols)) + 1; + end + end + %plot(pos(1,:), pos(2,:), '.') + %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, []); + pause + end +end + +eclass = bnet.equiv_class; +S=struct(bnet.CPD{eclass(Q2,2)}); + + + + + + + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/get_square_data.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/get_square_data.m new file mode 100644 index 00000000..9790221c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/get_square_data.m @@ -0,0 +1,70 @@ +% Let the user draw a square with the mouse, +% and then click on the corners to do a manual segmentation + +ss = 6; +Q1 = 1; Q2 = 2; Q3 = 3; obsvel = 6; +CLOCKWISE = 1; ANTICLOCK = 2; +LR = 1; UD = 2; RL = 3; DU = 4; + +% repeat this block manually incrementing the sequence number +% and setting ori. +% (since I don't know how to call getmouse as a call-return function). +seq = 4; +%ori = CLOCKWISE +ori = ANTICLOCK; +clear xpos ypos +getmouse +% end block + +% manual segmentation with the mouse +startseg(1) = 1; +for i=2:4 + fprintf('click on start of segment %d\n', i); + [x,y] = ginput(1); + plot(x,y,'ro') + d = dist2([xpos; ypos]', [x y]); + startseg(i) = argmin(d); +end + +% plot corners in green +%ti = first point in (i+1)st segment +t1 = startseg(1); t2 = startseg(2); t3 = startseg(3); t4 = startseg(4); +plot(xpos(t2), ypos(t2), 'g*') +plot(xpos(t3), ypos(t3), 'g*') +plot(xpos(t4), ypos(t4), 'g*') + + +xvel = xpos(2:end) - xpos(1:end-1); +yvel = ypos(2:end) - ypos(1:end-1); +speed = [xvel(:)'; yvel(:)']; +pos_data{seq} = [xpos(:)'; ypos(:)']; +vel_data{seq} = [xvel(:)'; yvel(:)']; +T = length(xvel); +Q1label{seq} = num2cell(repmat(ori, 1, T)); +Q2label{seq} = zeros(1, T); +if ori == CLOCKWISE + Q2label{seq}(t1:t2) = LR; + Q2label{seq}(t2+1:t3) = UD; + Q2label{seq}(t3+1:t4) = RL; + Q2label{seq}(t4+1:T) = DU; +else + Q2label{seq}(t1:t2) = RL; + Q2label{seq}(t2+1:t3) = UD; + Q2label{seq}(t3+1:t4) = LR; + Q2label{seq}(t4+1:T) = DU; +end + +% pos_data{seq}(:,t), vel_data{seq}(:,t) Q1label{seq}(t) Q2label{seq}(t) +save 'square4' pos_data vel_data Q1label Q2label + +nseq = 4; +cases = cell(1,nseq); +for seq=1:nseq + T = size(vel_data{seq},2); + ev = cell(ss,T); + ev(obsvel,:) = num2cell(vel_data{seq},1); + ev(Q1,:) = Q1label{seq}; + ev(Q2,:) = num2cell(Q2label{seq}); + cases{seq} = ev; +end +save 'square4_cases' cases diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/hhmm_inference.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/hhmm_inference.m new file mode 100644 index 00000000..c3bc8441 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/hhmm_inference.m @@ -0,0 +1,13 @@ +bnet = mk_square_hhmm(1, 1); + +engine = {}; +engine{end+1} = hmm_inf_engine(bnet); +engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet)); + +exact = 1:length(engine); +filter = 0; +single = 0; +maximize = 0; +T = 4; + +[err, inf_time, engine] = cmp_inference(bnet, engine, exact, T, filter, single, maximize); diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/is_F2_true_D3.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/is_F2_true_D3.m new file mode 100644 index 00000000..38d0b6e8 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/is_F2_true_D3.m @@ -0,0 +1,12 @@ +function stop = is_F2_true_D3(vals) +% function stop = is_F2_true_D3(vals) +% +% If vals(F2)=2 then level 2 has finished, so we return stop=1 +% to stop sample_dbn. Otherwise we return stop=0. +% We assume this is for a D=3 level HHMM. + +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +stop = 0; +if (iscell(vals) & vals{F2}==2) | (~iscell(vals) & vals(F2)==2) + stop = 1; +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_cts.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_cts.m new file mode 100644 index 00000000..77bdaec3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_cts.m @@ -0,0 +1,152 @@ +% Try to learn a 3 level HHMM similar to mk_square_hhmm +% from hand-drawn squares. + +% Because startprob should be shared for t=1:T, +% but in the DBN is shared for t=2:T, we train using a single long sequence. + +discrete_obs = 0; +supervised = 1; +obs_finalF2 = 0; +% It is not possible to observe F2 if we learn +% because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume +% the F nodes are always hidden (for speed). +% However, for generating, we might want to set the final F2=true +% to force all subroutines to finish. + +seed = 1; +rand('state', seed); +randn('state', seed); + +bnet = mk_square_hhmm(discrete_obs, 0); + +ss = 6; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; +Qsizes = [2 4 1]; + +if supervised + bnet.observed = [Q1 Q2 Onode]; +else + bnet.observed = [Onode]; +end + +if obs_finalF2 + engine = jtree_dbn_inf_engine(bnet); + % can't use ndx version because sometimes F2 is hidden, sometimes observed + error('can''t observe F when learning') +else + if supervised + engine = jtree_ndx_dbn_inf_engine(bnet); + else + engine = jtree_hmm_inf_engine(bnet); + end +end + +load 'square4_cases' % cases{seq}{i,t} for i=1:ss +%plot_square_hhmm(cases{1}) +%long_seq = cat(2, cases{:}); +train_cases = cases(1:2); +long_seq = cat(2, train_cases{:}); +if ~supervised + T = size(long_seq,2); + for t=1:T + long_seq{Q1,t} = []; + long_seq{Q2,t} = []; + end +end +[bnet2, LL, engine2] = learn_params_dbn_em(engine, {long_seq}, 'max_iter', 2); + +eclass = bnet2.equiv_class; +CPDO=struct(bnet2.CPD{eclass(Onode,1)}); +mu = CPDO.mean; +Sigma = CPDO.cov; +CPDO_full = CPDO; + +% force diagonal covs after training +for k=1:size(Sigma,3) + Sigma(:,:,k) = diag(diag(Sigma(:,:,k))); +end +bnet2.CPD{6} = set_fields(bnet.CPD{6}, 'cov', Sigma); + +if 0 + % visualize each model by concatenating means for each model for nsteps in a row + nsteps = 5; + ev = cell(ss, nsteps*prod(Qsizes(2:3))); + t = 1; + for q2=1:Qsizes(2) + for q3=1:Qsizes(3) + for i=1:nsteps + ev{Onode,t} = mu(:,q2,q3); + ev{Q2,t} = q2; + t = t + 1; + end + end + end + plot_square_hhmm(ev) +end + +% bnet3 is the same as the learned model, except we will use it in testing mode +if supervised + bnet3 = bnet2; + bnet3.observed = [Onode]; + engine3 = hmm_inf_engine(bnet3); + %engine3 = jtree_ndx_dbn_inf_engine(bnet3); +else + bnet3 = bnet2; + engine3 = engine2; +end + +if 0 + % segment whole sequence + mpe = calc_mpe_dbn(engine3, long_seq); + pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []); +end + +% segment each sequence +test_cases = cases(3:4); +for i=1:2 + ev = test_cases{i}; + T = size(ev, 2); + for t=1:T + ev{Q1,t} = []; + ev{Q2,t} = []; + end + %mpe = calc_mpe_dbn(engine3, ev); + mpe = find_mpe(engine3, ev) + subplot(1,2,i) + plot_square_hhmm(mpe) + %pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, []); + q1s = cell2num(mpe(Q1,:)); + h = hist(q1s, 1:Qsizes(1)); + map_q1 = argmax(h); + str = sprintf('test seq %d is of type %d\n', i, map_q1); + title(str) +end + + +if 0 +% Estimate gotten by couting transitions in the labelled data +% Note that a self transition shouldnt count if F2=off. +Q2ev = cell2num(ev(Q2,:)); +Q2a = Q2ev(1:end-1); +Q2b = Q2ev(2:end); +counts = compute_counts([Q2a; Q2b], [4 4]); +end + +eclass = bnet2.equiv_class; +CPDQ1=struct(bnet2.CPD{eclass(Q1,2)}); +CPDQ2=struct(bnet2.CPD{eclass(Q2,2)}); +CPDQ3=struct(bnet2.CPD{eclass(Q3,2)}); +CPDF2=struct(bnet2.CPD{eclass(F2,1)}); +CPDF3=struct(bnet2.CPD{eclass(F3,1)}); + + +A=add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2)); +squeeze(A(:,1,:)); +CPDQ2.startprob; + +if 0 +S=struct(CPDF2.sub_CPD_term); +S.nsamples +reshape(S.counts, [2 4 2]) +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_discrete.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_discrete.m new file mode 100644 index 00000000..3110ae5e --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/learn_square_hhmm_discrete.m @@ -0,0 +1,171 @@ +% Try to learn a 3 level HHMM similar to mk_square_hhmm +% from synthetic discrete sequences + + +discrete_obs = 1; +supervised = 0; +obs_finalF2 = 0; + +seed = 1; +rand('state', seed); +randn('state', seed); + +bnet_init = mk_square_hhmm(discrete_obs, 0); + +ss = 6; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; + +if supervised + bnet_init.observed = [Q1 Q2 Onode]; +else + bnet_init.observed = [Onode]; +end + +if obs_finalF2 + engine_init = jtree_dbn_inf_engine(bnet_init); + % can't use ndx version because sometimes F2 is hidden, sometimes observed + error('can''t observe F when learning') + % It is not possible to observe F2 if we learn + % because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume + % the F nodes are always hidden (for speed). + % However, for generating, we might want to set the final F2=true + % to force all subroutines to finish. +else + if supervised + engine_init = jtree_ndx_dbn_inf_engine(bnet_init); + else + engine_init = hmm_inf_engine(bnet_init); + end +end + +% generate some synthetic data (easier to debug) +chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; +L=find(chars=='L'); l=find(chars=='l'); +U=find(chars=='U'); u=find(chars=='u'); +R=find(chars=='R'); r=find(chars=='r'); +D=find(chars=='D'); d=find(chars=='d'); + +cases = {}; + +T = 8; +ev = cell(ss, T); +ev(Onode,:) = num2cell([L l U u R r D d]); +if supervised + ev(Q1,:) = num2cell(1*ones(1,T)); + ev(Q2,:) = num2cell( [1 1 2 2 3 3 4 4]); +end +cases{1} = ev; +cases{3} = ev; + +T = 8; +ev = cell(ss, T); +%we start with R then r, even though we are running the model 'backwards'! +ev(Onode,:) = num2cell([R r U u L l D d]); + +if supervised + ev(Q1,:) = num2cell(2*ones(1,T)); + ev(Q2,:) = num2cell( [3 3 2 2 1 1 4 4]); +end + +cases{2} = ev; +cases{4} = ev; + +if obs_finalF2 + for i=1:length(cases) + T = size(cases{i},2); + cases{i}(F2,T)={2}; % force F2 to be finished at end of seq + end +end + + +% startprob should be shared for t=1:T, +% but in the DBN it is shared for t=2:T, +% so we train using a single long sequence. +long_seq = cat(2, cases{:}); +[bnet_learned, LL, engine_learned] = ... + learn_params_dbn_em(engine_init, {long_seq}, 'max_iter', 200); + +% figure out which subsequence each model is responsible for +mpe = calc_mpe_dbn(engine_learned, long_seq); +pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, chars); + + +% The "true" segmentation of the training sequence is +% Q1: 1 2 +% O: L l U u R r D d | R r U u L l D d | etc. +% +% When we learn in a supervised fashion, we recover the "truth". + +% When we learn in an unsupervised fashion with seed=1, we get +% Q1: 2 1 +% O: L l U u R r D d R r | U u L l D d | etc. +% +% This means for model 1: +% starts in state 2 +% transitions 2->1, 1->4, 4->e, 3->2 +% +% For model 2, +% starts in state 1 +% transitions 1->2, 2->3, 3->4 or e, 4->3 + +% examine the params +eclass = bnet_learned.equiv_class; +CPDQ1=struct(bnet_learned.CPD{eclass(Q1,2)}); +CPDQ2=struct(bnet_learned.CPD{eclass(Q2,2)}); +CPDQ3=struct(bnet_learned.CPD{eclass(Q3,2)}); +CPDF2=struct(bnet_learned.CPD{eclass(F2,1)}); +CPDF3=struct(bnet_learned.CPD{eclass(F3,1)}); +CPDO=struct(bnet_learned.CPD{eclass(Onode,1)}); + +A_learned =add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2)); +squeeze(A_learned(:,1,:)) +squeeze(A_learned(:,2,:)) + + +% Does the "true" model have higher likelihood than the learned one? +% i.e., Does the unsupervised method learn the wrong model because +% we have the wrong cost fn, or because of local minima? + +bnet_true = mk_square_hhmm(discrete_obs,1); + +% examine the params +eclass = bnet_learned.equiv_class; +CPDQ1_true=struct(bnet_true.CPD{eclass(Q1,2)}); +CPDQ2_true=struct(bnet_true.CPD{eclass(Q2,2)}); +CPDQ3_true=struct(bnet_true.CPD{eclass(Q3,2)}); +CPDF2_true=struct(bnet_true.CPD{eclass(F2,1)}); +CPDF3_true=struct(bnet_true.CPD{eclass(F3,1)}); + +A_true =add_hhmm_end_state(CPDQ2_true.transprob, CPDF2_true.termprob(:,:,2)); +squeeze(A_true(:,1,:)) + + +if supervised + engine_true = jtree_ndx_dbn_inf_engine(bnet_true); +else + engine_true = hmm_inf_engine(bnet_true); +end + +%[engine_learned, ll_learned] = enter_evidence(engine_learned, long_seq); +%[engine_true, ll_true] = enter_evidence(engine_true, long_seq); +[engine_learned, ll_learned] = enter_evidence(engine_learned, cases{2}); +[engine_true, ll_true] = enter_evidence(engine_true, cases{2}); +ll_learned +ll_true + + +% remove concatentation artefacts +ll_learned = 0; +ll_true = 0; +for m=1:length(cases) + [engine_learned, ll_learned_tmp] = enter_evidence(engine_learned, cases{m}); + [engine_true, ll_true_tmp] = enter_evidence(engine_true, cases{m}); + ll_learned = ll_learned + ll_learned_tmp; + ll_true = ll_true + ll_true_tmp; +end +ll_learned +ll_true + +% In both cases, ll_learned >> ll_true +% which shows we are using the wrong cost function! diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/mk_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/mk_square_hhmm.m new file mode 100644 index 00000000..41cfc539 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/mk_square_hhmm.m @@ -0,0 +1,180 @@ +function bnet = mk_square_hhmm(discrete_obs, true_params, topright) + +% Make a 3 level HHMM described by the following grammar +% +% Square -> CLK | CCK % clockwise or counterclockwise +% CLK -> LR UD RL DU start on top left (1 2 3 4) +% CCK -> RL UD LR DU if start at top right (3 2 1 4) +% CCK -> UD LR DU RL if start at top left (2 1 4 3) +% +% LR = left-right, UD = up-down, RL = right-left, DU = down-up +% LR, UD, RL, DU are sub HMMs. +% +% For discrete observations, the subHMMs are 2-state left-right. +% LR emits L then l, etc. +% +% For cts observations, the subHMMs are 1 state. +% LR emits a vector in the -> direction, with a little noise. +% Since there is no constraint that we remain in the LR state as long as the RL state, +% the sides of the square might have different lengths, +% so the result is not really a square! +% +% If true_params = 0, we use random parameters at the top 2 levels +% (ready for learning). At the bottom level, we use noisy versions +% of the "true" observations. +% +% If topright=1, counter-clockwise starts at top right, not top left +% This example was inspired by Ivanov and Bobick. + +if nargin < 3, topright = 1; end + +if 1 % discrete_obs + Qsizes = [2 4 2]; +else + Qsizes = [2 4 1]; +end + +D = 3; +Qnodes = 1:D; +startprob = cell(1,D); +transprob = cell(1,D); +termprob = cell(1,D); + +% LEVEL 1 + +startprob{1} = 'unif'; +transprob{1} = 'unif'; + +% LEVEL 2 + +if true_params + startprob{2} = zeros(2, 4); + startprob{2}(1, :) = [1 0 0 0]; + if topright + startprob{2}(2, :) = [0 0 1 0]; + else + startprob{2}(2, :) = [0 1 0 0]; + end + + transprob{2} = zeros(4, 2, 4); + + transprob{2}(:,1,:) = [0 1 0 0 + 0 0 1 0 + 0 0 0 1 + 0 0 0 1]; % 4->e + if topright + transprob{2}(:,2,:) = [0 0 0 1 + 1 0 0 0 + 0 1 0 0 + 0 0 0 1]; % 4->e + else + transprob{2}(:,2,:) = [0 0 0 1 + 1 0 0 0 + 0 0 1 0 % 3->e + 0 0 1 0]; + end + + %termprob{2} = 'rightstop'; + termprob{2} = zeros(2,4); + pfin = 0.8; + termprob{2}(1,:) = [0 0 0 pfin]; % finish in state 4 (DU) + if topright + termprob{2}(2,:) = [0 0 0 pfin]; + else + termprob{2}(2,:) = [0 0 pfin 0]; % finish in state 3 (RL) + end +else + % In the unsupervised case, it is essential that we break symmetry + % in the initial param estimates. + %startprob{2} = 'unif'; + %transprob{2} = 'unif'; + %termprob{2} = 'unif'; + startprob{2} = 'rnd'; + transprob{2} = 'rnd'; + termprob{2} = 'rnd'; +end + +% LEVEL 3 + +if 1 | true_params + startprob{3} = 'leftstart'; + transprob{3} = 'leftright'; + termprob{3} = 'rightstop'; +else + % If we want to be able to run a base-level model backwards... + startprob{3} = 'rnd'; + transprob{3} = 'rnd'; + termprob{3} = 'rnd'; +end + + +% OBS LEVEl + +if discrete_obs + % Initialise observations of lowest level primitives in a way which we can interpret + chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; + L=find(chars=='L'); l=find(chars=='l'); + U=find(chars=='U'); u=find(chars=='u'); + R=find(chars=='R'); r=find(chars=='r'); + D=find(chars=='D'); d=find(chars=='d'); + Osize = length(chars); + + if true_params + p = 1; % makes each state fully observed + else + p = 0.9; + end + + obsprob = (1-p)*ones([4 2 Osize]); + % Q2 Q3 O + obsprob(1, 1, L) = p; + obsprob(1, 2, l) = p; + obsprob(2, 1, U) = p; + obsprob(2, 2, u) = p; + obsprob(3, 1, R) = p; + obsprob(3, 2, r) = p; + obsprob(4, 1, D) = p; + obsprob(4, 2, d) = p; + obsprob = mk_stochastic(obsprob); + Oargs = {'CPT', obsprob}; +else + % Initialise means of lowest level primitives in a way which we can interpret + % These means are little vectors in the east, south, west, north directions. + % (left-right=east, up-down=south, right-left=west, down-up=north) + Osize = 2; + mu = zeros(2, Qsizes(2), Qsizes(3)); + scale = 3; + if true_params + noise = 0; + else + noise = 0.5*scale; + end + for q3=1:Qsizes(3) + mu(:, 1, q3) = scale*[1;0] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 2, q3) = scale*[0;-1] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 3, q3) = scale*[-1;0] + noise*rand(2,1); + end + for q3=1:Qsizes(3) + mu(:, 4, q3) = scale*[0;1] + noise*rand(2,1); + end + Sigma = repmat(reshape(scale*eye(2), [2 2 1 1 ]), [1 1 Qsizes(2) Qsizes(3)]); + Oargs = {'mean', mu, 'cov', Sigma, 'cov_type', 'diag'}; +end + +if discrete_obs + selfprob = 0.5; +else + selfprob = 0.95; + % If less than this, it won't look like a square + % because it doesn't spend enough time in each state + % Unfortunately, the variance on durations (lengths of each side) + % is very large +end +bnet = mk_hhmm('Qsizes', Qsizes, 'Osize', Osize', 'discrete_obs', discrete_obs, ... + 'Oargs', Oargs, 'Ops', Qnodes(2:3), 'selfprob', selfprob, ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/plot_square_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/plot_square_hhmm.m new file mode 100644 index 00000000..e61e5669 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/plot_square_hhmm.m @@ -0,0 +1,27 @@ +function plot_square_hhmm(ev) +% Plot the square shape implicit in the evidence. +% ev{i,t} is the value of node i in slice t. +% The observed node contains a velocity (delta increment), which is converted +% into a position. +% The Q2 node specifies which model is used, and hence which color +% to use: 1=red, 2=green, 3=blue, 4=black. + +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; + +delta = cell2num(ev(Onode,:)); % delta(:,t) +Q2label = cell2num(ev(Q2,:)); + +T = size(delta, 2); +pos = zeros(2,T+1); +hold on +cols = {'r', 'g', 'b', 'k'}; +for t=2:T+1 + pos(:,t) = pos(:,t-1) + delta(:,t-1); + plot(pos(1,t), pos(2,t), sprintf('%c.', cols{Q2label(t-1)})); + if (t==2) + text(pos(1,t-1),pos(2,t-1),sprintf('%d',t)) + elseif (mod(t,20)==0) + text(pos(1,t),pos(2,t),sprintf('%d',t)) + end +end + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_cts.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_cts.m new file mode 100644 index 00000000..3ab2abb5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_cts.m @@ -0,0 +1,20 @@ +% Generate samples from the HHMM with the true params. + +seed = 1; +rand('state', seed); +randn('state', seed); + +discrete_obs = 0; + +bnet = mk_square_hhmm(discrete_obs, 1); +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; + +for seqi=1:1 + evidence = sample_dbn(bnet, 'stop_test', 'is_F2_true_D3'); + clf + plot_square_hhmm(evidence); + %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, []); + fprintf('sequence %d has length %d; press key to continue\n', seqi, size(evidence,2)) + pause +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_discrete.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_discrete.m new file mode 100644 index 00000000..20279899 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/sample_square_hhmm_discrete.m @@ -0,0 +1,20 @@ +% Generate samples from the HHMM with the true params. + +seed = 0; +rand('state', seed); +randn('state', seed); + +discrete_obs = 1; + +bnet = mk_square_hhmm(discrete_obs, 1); + +Tmax = 30; +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; +chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd']; + +for seqi=1:3 + evidence = cell2num(sample_dbn(bnet, 'stop_test', 'is_F2_true_D3')); + T = size(evidence, 2) + pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars); +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4.mat b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4.mat new file mode 100644 index 00000000..cda0585b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4.mat Binary files differdiff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4_cases.mat b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4_cases.mat new file mode 100644 index 00000000..788c3239 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/square4_cases.mat Binary files differdiff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.m new file mode 100644 index 00000000..e983af14 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.m @@ -0,0 +1,1310 @@ +function fig = test_square_fig() +% This is the machine-generated representation of a Handle Graphics object +% and its children. Note that handle values may change when these objects +% are re-created. This may cause problems with any callbacks written to +% depend on the value of the handle at the time the object was saved. +% +% To reopen this object, just type the name of the M-file at the MATLAB +% prompt. The M-file and its associated MAT-file must be on your path. + +load test_square_fig + +h0 = figure('Color',[0.8 0.8 0.8], ... + 'Colormap',mat0, ... + 'PointerShapeCData',mat1, ... + 'Position',[540 374 476 292]); +h1 = axes('Parent',h0, ... + 'CameraUpVector',[0 1 0], ... + 'Color',[1 1 1], ... + 'ColorOrder',mat2, ... + 'NextPlot','add', ... + 'Position',[0.13 0.11 0.3270231213872832 0.8149999999999998], ... + 'XColor',[0 0 0], ... + 'XLim',[-10 50], ... + 'XLimMode','manual', ... + 'YColor',[0 0 0], ... + 'YLim',[-60 10], ... + 'YLimMode','manual', ... + 'ZColor',[0 0 0]); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',0.4608294930875587, ... + 'YData',0.2923976608187218); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'String','2'); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',1.152073732718893, ... + 'YData',0.2923976608187218); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',2.995391705069125, ... + 'YData',0.8771929824561511); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',3.686635944700463, ... + 'YData',0.8771929824561511); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',6.451612903225808, ... + 'YData',0.8771929824561511); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',9.677419354838712, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',10.36866359447005, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',15.43778801843318, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',17.51152073732719, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',19.81566820276498, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',20.50691244239631, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',23.73271889400922, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',25.57603686635945, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',29.95391705069125, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',31.79723502304147, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',35.02304147465438, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',35.71428571428572, ... + 'YData',2.046783625730996); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',38.47926267281106, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',40.3225806451613, ... + 'YData',1.461988304093566); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[40.3225806451613 1.461988304093566 0], ... + 'String','20'); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',42.62672811059908, ... + 'YData',mat3); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',43.31797235023042, ... + 'YData',mat4); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',43.31797235023042, ... + 'YData',0.8771929824561511); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',43.54838709677419, ... + 'YData',0); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',43.77880184331798, ... + 'YData',-0.5847953216374293); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',44.47004608294931, ... + 'YData',-2.339181286549703); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',44.93087557603687, ... + 'YData',-4.385964912280699); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',46.7741935483871, ... + 'YData',-9.064327485380119); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.00460829493088, ... + 'YData',-10.81871345029239); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.69585253456221, ... + 'YData',mat5); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.69585253456221, ... + 'YData',-15.20467836257309); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-19.00584795321637); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-19.88304093567251); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-22.51461988304093); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-23.09941520467836); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-26.02339181286549); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-26.31578947368421); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-27.77777777777777); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-28.3625730994152); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.23502304147466, ... + 'YData',-30.99415204678362); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[47.23502304147466 -30.99415204678362 0], ... + 'String','40'); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.46543778801843, ... + 'YData',-31.57894736842105); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.46543778801843, ... + 'YData',-33.62573099415204); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.46543778801843, ... + 'YData',-34.50292397660818); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.46543778801843, ... + 'YData',-37.42690058479531); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.46543778801843, ... + 'YData',-38.01169590643274); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.00460829493088, ... + 'YData',-42.39766081871344); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.00460829493088, ... + 'YData',-42.98245614035087); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',47.00460829493088, ... + 'YData',-46.49122807017543); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',46.7741935483871, ... + 'YData',-46.78362573099415); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',46.54377880184332, ... + 'YData',-49.41520467836257); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',46.54377880184332, ... + 'YData',-49.70760233918128); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',45.85253456221199, ... + 'YData',-51.46198830409356); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',45.85253456221199, ... + 'YData',-51.75438596491227); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',44.93087557603687, ... + 'YData',-53.21637426900584); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',44.70046082949308, ... + 'YData',-53.21637426900584); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',44.00921658986175, ... + 'YData',-54.09356725146198); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',43.77880184331798, ... + 'YData',-54.38596491228069); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',41.93548387096774, ... + 'YData',-54.97076023391811); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',41.47465437788019, ... + 'YData',-55.26315789473683); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',39.1705069124424, ... + 'YData',-55.55555555555554); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[39.1705069124424 -55.55555555555554 0], ... + 'String','60'); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',38.94009216589862, ... + 'YData',-55.84795321637426); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',36.63594470046083, ... + 'YData',-55.55555555555554); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',36.17511520737327, ... + 'YData',-55.55555555555554); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',32.94930875576037, ... + 'YData',-54.97076023391811); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',32.48847926267281, ... + 'YData',-54.97076023391811); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',28.11059907834102, ... + 'YData',-53.80116959064326); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',27.64976958525346, ... + 'YData',-53.50877192982455); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',23.963133640553, ... + 'YData',-53.50877192982455); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',23.27188940092166, ... + 'YData',-53.50877192982455); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',19.5852534562212, ... + 'YData',-54.97076023391811); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',19.12442396313364, ... + 'YData',-54.97076023391811); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',mat6, ... + 'YData',-56.14035087719297); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',mat7, ... + 'YData',-56.14035087719297); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',9.907834101382491, ... + 'YData',-57.30994152046782); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',9.447004608294932, ... + 'YData',-57.30994152046782); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',6.221198156682029, ... + 'YData',-57.30994152046782); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',4.838709677419356, ... + 'YData',-56.7251461988304); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',2.764976958525345, ... + 'YData',-56.14035087719297); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',2.534562211981569, ... + 'YData',-56.14035087719297); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',0.9216589861751174, ... + 'YData',-53.80116959064327); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[0.9216589861751174 -53.80116959064327 0], ... + 'String','80'); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',0.6912442396313381, ... + 'YData',-53.21637426900584); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.152073732718893, ... + 'YData',-48.24561403508771); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.152073732718893, ... + 'YData',-47.953216374269); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.843317972350228, ... + 'YData',-44.73684210526315); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.843317972350228, ... + 'YData',-44.44444444444444); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-2.304147465437787, ... + 'YData',-39.76608187134502); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-2.764976958525345, ... + 'YData',-38.01169590643274); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-3.225806451612904, ... + 'YData',-30.99415204678362); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-3.225806451612904, ... + 'YData',-29.82456140350877); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-3.225806451612904, ... + 'YData',-24.85380116959064); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-3.225806451612904, ... + 'YData',-24.26900584795321); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-2.534562211981566, ... + 'YData',-17.5438596491228); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-2.304147465437787, ... + 'YData',-16.95906432748537); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.612903225806452, ... + 'YData',-11.98830409356725); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.612903225806452, ... + 'YData',-11.40350877192982); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',mat8, ... + 'YData',-8.47953216374269); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',mat9, ... + 'YData',-8.187134502923968); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.152073732718893, ... + 'YData',-5.263157894736835); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-1.152073732718893, ... + 'YData',-4.970760233918128); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.9216589861751139, ... + 'YData',-2.923976608187132); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[-0.9216589861751139 -2.923976608187132 0], ... + 'String','100'); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.9216589861751139, ... + 'YData',-2.631578947368411); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.6912442396313345, ... + 'YData',mat10); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.6912442396313345, ... + 'YData',-0.8771929824561369); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.6912442396313345, ... + 'YData',-0.5847953216374293); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'HandleVisibility','off', ... + 'HorizontalAlignment','center', ... + 'Position',[19.80645161290322 12.0675105485232 17.32050807568877], ... + 'VerticalAlignment','bottom'); +set(get(h2,'Parent'),'Title',h2); +h1 = axes('Parent',h0, ... + 'CameraUpVector',[0 1 0], ... + 'Color',[1 1 1], ... + 'ColorOrder',mat11, ... + 'NextPlot','add', ... + 'Position',[0.5779768786127169 0.11 0.3270231213872832 0.8149999999999998], ... + 'XColor',[0 0 0], ... + 'YColor',[0 0 0], ... + 'ZColor',[0 0 0]); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.4608294930875587, ... + 'YData',0); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'String','2'); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-2.764976958525345, ... + 'YData',-0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-3.456221198156683, ... + 'YData',-0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-7.834101382488477, ... + 'YData',-0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-11.52073732718894, ... + 'YData',-0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',mat12, ... + 'YData',-0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-19.35483870967742, ... + 'YData',0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-23.50230414746544, ... + 'YData',0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-24.88479262672811, ... + 'YData',0.8771929824561369); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-28.11059907834102, ... + 'YData',0.8771929824561369); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-29.49308755760369, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-31.10599078341014, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-32.02764976958525, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-33.17972350230414, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-33.6405529953917, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 1], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-34.7926267281106, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-35.02304147465438, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-35.48387096774194, ... + 'YData',1.461988304093566); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-35.71428571428572, ... + 'YData',1.461988304093566); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[-35.71428571428572 1.461988304093566 0], ... + 'String','20'); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.17511520737327, ... + 'YData',mat13); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',0.8771929824561369); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.63594470046083, ... + 'YData',0.2923976608187076); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.63594470046083, ... + 'YData',0); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.63594470046083, ... + 'YData',mat14); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',-2.339181286549703); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',-2.631578947368425); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',-4.67836257309942); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',-5.555555555555557); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',-8.187134502923982); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.40552995391705, ... + 'YData',-8.771929824561397); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.86635944700461, ... + 'YData',-13.15789473684211); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',mat15); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-16.08187134502924); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-17.54385964912281); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-18.12865497076023); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-19.88304093567251); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-20.17543859649123); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-21.92982456140351); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-22.22222222222222); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[-37.09677419354839 -22.22222222222222 0], ... + 'String','40'); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-23.09941520467836); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.09677419354839, ... + 'YData',-23.39181286549707); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-25.14619883040935); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-25.43859649122807); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-28.3625730994152); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.55760368663595, ... + 'YData',-28.94736842105263); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-31.87134502923976); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-32.16374269005848); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-34.7953216374269); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-35.38011695906432); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-38.88888888888889); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-39.76608187134503); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-43.27485380116958); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-43.5672514619883); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.78801843317973, ... + 'YData',-44.44444444444444); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.55760368663595, ... + 'YData',-45.32163742690058); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.55760368663595, ... + 'YData',-45.61403508771929); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-47.36842105263158); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-37.32718894009217, ... + 'YData',-47.95321637426901); +h2 = line('Parent',h1, ... + 'Color',[0 1 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.86635944700461, ... + 'YData',-49.70760233918129); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[-36.86635944700461 -49.70760233918129 0], ... + 'String','60'); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-36.86635944700461, ... + 'YData',-50); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-35.71428571428572, ... + 'YData',-50.29239766081872); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-35.25345622119816, ... + 'YData',-50.29239766081872); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-32.02764976958527, ... + 'YData',-50.29239766081872); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-31.33640552995393, ... + 'YData',-50.29239766081872); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-27.88018433179725, ... + 'YData',-50.58479532163743); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-27.41935483870969, ... + 'YData',-50.58479532163743); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-18.20276497695854, ... + 'YData',-50.58479532163743); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-16.82027649769586, ... + 'YData',-51.16959064327486); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-12.21198156682029, ... + 'YData',-50.58479532163743); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-11.52073732718895, ... + 'YData',-50.58479532163743); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-6.912442396313377, ... + 'YData',-51.16959064327486); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-5.069124423963142, ... + 'YData',-51.75438596491229); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',mat16, ... + 'YData',-52.046783625731); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',-0.9216589861751281, ... + 'YData',-52.33918128654972); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',0.2304147465437687, ... + 'YData',-52.33918128654972); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',0.4608294930875481, ... + 'YData',-52.33918128654972); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',2.304147465437776, ... + 'YData',-52.63157894736843); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',2.534562211981548, ... + 'YData',-52.63157894736843); +h2 = line('Parent',h1, ... + 'Color',[1 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',3.917050691244224, ... + 'YData',-52.63157894736843); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[3.917050691244224 -52.63157894736843 0], ... + 'String','80'); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',4.147465437788011, ... + 'YData',-52.63157894736843); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',4.147465437788011, ... + 'YData',-52.33918128654972); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.529953917050673, ... + 'YData',-46.19883040935674); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.990783410138231, ... + 'YData',-44.44444444444446); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',7.834101382488466, ... + 'YData',-28.0701754385965); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',8.294930875576025, ... + 'YData',-22.80701754385966); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',8.755760368663584, ... + 'YData',mat17); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',8.525345622119797, ... + 'YData',-14.9122807017544); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',7.373271889400908, ... + 'YData',-10.23391812865498); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',7.142857142857135, ... + 'YData',-9.94152046783627); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',6.221198156682018, ... + 'YData',-7.602339181286567); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',6.221198156682018, ... + 'YData',-7.309941520467845); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.760368663594459, ... + 'YData',-5.555555555555571); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.529953917050673, ... + 'YData',-5.555555555555571); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.2995391705069, ... + 'YData',-4.093567251462005); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.069124423963128, ... + 'YData',-2.631578947368439); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.069124423963128, ... + 'YData',-2.339181286549717); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.069124423963128, ... + 'YData',mat18); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',5.069124423963128, ... + 'YData',-1.169590643274873); +h2 = line('Parent',h1, ... + 'Color',[0 0 0], ... + 'LineStyle','none', ... + 'Marker','.', ... + 'XData',4.838709677419342, ... + 'YData',-0.2923976608187218); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'Position',[4.838709677419342 -0.2923976608187218 0], ... + 'String','100'); +h2 = text('Parent',h1, ... + 'Color',[0 0 0], ... + 'HandleVisibility','off', ... + 'HorizontalAlignment','center', ... + 'Position',[-10.38961038961038 12.0675105485232 17.32050807568877], ... + 'VerticalAlignment','bottom'); +set(get(h2,'Parent'),'Title',h2); +if nargout > 0, fig = h0; end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.mat b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.mat new file mode 100644 index 00000000..5b2b5f53 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Square/test_square_fig.mat Binary files differdiff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/abcd_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/abcd_hhmm.m new file mode 100644 index 00000000..13e6d39f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/abcd_hhmm.m @@ -0,0 +1,97 @@ +% Make the HHMM in Figure 1 of the NIPS'01 paper + +Qsize = [2 3 2]; +Qnodes = 1:3; +D = 3; +transprob = cell(1,D); +termprob = cell(1,D); +startprob = cell(1,D); +clear A; + +% transprob{d}(i,k,j), transprob{1}(i,j) +% termprob{d}(k,j), termprob{1}(1,j) +% startprob{d}(k,j), startprob{1}(1,j) + + +% LEVEL 1 + +% 1 2 e +A{1} = [0 0 1; + 0 0 1]; +[transprob{1}, termprob{1}] = remove_hhmm_end_state(A{1}); +startprob{1} = [0.5 0.5]; + +% LEVEL 2 +A{2} = zeros(Qsize(2), Qsize(1), Qsize(2)+1); + +% 1 2 3 e +A{2}(:,1,:) = [0 1 0 0 % Q1=1 => model below state 0 + 0 0 1 0 + 0 0 0 1]; + +% 1 2 3 e +A{2}(:,2,:) = [0 1 0 0 % Q1=2 => model below state 1 + 0 0 1 0 + 0 0 0 1]; + +[transprob{2}, termprob{2}] = remove_hhmm_end_state(A{2}); + +% always enter level 2 in state 1 +startprob{2} = [1 0 0 + 1 0 0]; + +% LEVEL 3 + +A{3} = zeros([Qsize(3) Qsize(2) Qsize(3)+1]); +endstate = Qsize(3)+1; +% Qt-1(3) Qt(2) Qt(3) +% 1 2 e +A{3}(1, 1, endstate) = 1.0; % Q2=1 => model below state 2/5 +A{3}(:, 2, :) = [0.0 1.0 0.0 % Q2=2 => model below state 3/6 + 0.5 0.0 0.5]; +A{3}(1, 3, endstate) = 1.0; % Q2=3 => model below state 4/7 + +[transprob{3}, termprob{3}] = remove_hhmm_end_state(A{3}); + +startprob{3} = 'leftstart'; + + + +% OBS LEVEl + +chars = ['a', 'b', 'c', 'd', 'x', 'y']; +Osize = length(chars); + +obsprob = zeros([Qsize Osize]); +% 1 2 3 O +obsprob(1,1,1,find(chars == 'a')) = 1.0; + +obsprob(1,2,1,find(chars == 'x')) = 1.0; +obsprob(1,2,2,find(chars == 'y')) = 1.0; + +obsprob(1,3,1,find(chars == 'b')) = 1.0; + +obsprob(2,1,1,find(chars == 'c')) = 1.0; + +obsprob(2,2,1,find(chars == 'x')) = 1.0; +obsprob(2,2,2,find(chars == 'y')) = 1.0; + +obsprob(2,3,1,find(chars == 'd')) = 1.0; + +Oargs = {'CPT', obsprob}; + +bnet = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ... + 'Oargs', Oargs, 'Ops', Qnodes(1:3), ... + 'startprob', startprob, 'transprob', transprob, 'termprob', termprob); + + +Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6; +Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3]; + +for seqi=1:3 + evidence = sample_dbn(bnet, 'stop_test', 'is_F2_true_D3'); + ev = cell2num(evidence); + chars(ev(end,:)) + %T = size(evidence, 2) + %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars); +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/add_hhmm_end_state.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/add_hhmm_end_state.m new file mode 100644 index 00000000..84c3c653 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/add_hhmm_end_state.m @@ -0,0 +1,34 @@ +function A = add_hhmm_end_state(transprob, termprob) +% ADD_HMM_END_STATE Combine trans and term probs into transmat for automaton with an end state +% function A = add_hhmm_end_state(transprob, termprob) +% +% A(i,k,j) = Pr( i->j | Qps=k), where i in 1:Q, j in 1:(Q+1), and Q+1 is the end state +% This implements the equation in sec 4.6 of my tech report, where +% transprob(i,k,j) = \tilde{A}_k(i,j), termprob(k,j) = \tau_k(j) +% +% For the top level, the k index is missing. + +Q = size(transprob,1); +toplevel = (ndims(transprob)==2); +if toplevel + Qk = 1; + transprob = reshape(transprob, [Q 1 Q]); + termprob = reshape(termprob, [1 Q]); +else + Qk = size(transprob, 2); +end + +A = zeros(Q, Qk, Q+1); +A(:,:,Q+1) = termprob'; + +for k=1:Qk + for i=1:Q + for j=1:Q + A(i,k,j) = transprob(i,k,j) * (1-termprob(k,i)); + end + end +end + +if toplevel + A = squeeze(A); +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/hhmm_jtree_clqs.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/hhmm_jtree_clqs.m new file mode 100644 index 00000000..4192af11 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/hhmm_jtree_clqs.m @@ -0,0 +1,143 @@ +% Find out how big the cliques are in an HHMM as a function of depth +% (This is how we get the complexity bound of O(D K^{1.5D}).) + +if 0 +Qsize = []; +Fsize = []; +Nclqs = []; +end + +ds = 1:15; + +for d = ds + allQ = 1; + [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(d, allQ); + + N = length(intra); + ns = 2*ones(1,N); + + bnet = mk_dbn(intra, inter, ns); + for i=1:N + bnet.CPD{i} = tabular_CPD(bnet, i); + end + + if 0 + T = 5; + dag = unroll_dbn_topology(intra, inter, T); + engine = jtree_unrolled_dbn_inf_engine(bnet, T, 'constrained', 1); + S = struct(engine); + S1 = struct(S.sub_engine); + end + + engine = jtree_dbn_inf_engine(bnet); + S = struct(engine); + J = S.jtree_struct; + + ss = 2*d+1; + Qnodes2 = Qnodes + ss; + QQnodes = [Qnodes Qnodes2]; + + % find out how many Q nodes in each clique, and how many F nodes + C = length(J.cliques); + Nclqs(d) = 0; + for c=1:C + Qsize(c,d) = length(myintersect(J.cliques{c}, QQnodes)); + Fsize(c,d) = length(myintersect(J.cliques{c}, Fnodes)); + if length(J.cliques{c}) > 1 % exclude observed leaves + Nclqs(d) = Nclqs(d) + 1; + end + end + %pred_max_Qsize(d) = ceil(d+(d+1)/2); + pred_max_Qsize(d) = ceil(1.5*d); + + fprintf('d=%d\n', d); + %fprintf('D=%d, max F = %d. max Q = %d, pred max Q = %d\n', ... + % D, max(Fsize), max(Qsize), ceil(D+(D+1)/2)); + + %histc(Qsize,1:max(Qsize)) % how many of each size? +end % next d + + +Q = 2; +pred_mass = ds.*(Q.^ds) + Q.^(ceil(1.5 * ds)) +pred_mass2 = Q.^(ceil(1.5 * ds)) + +for d=ds + mass(d) = 0; + for c=1:C + mass(d) = mass(d) + Q^Qsize(c,d); + end +end + + +if 0 +%plot(ds, max(Qsize), 'o-', ds, pred_max_Qsize, '*--'); +%plot(ds, max(Qsize), 'o-', ds, 1.5*ds, '*--'); +%plot(ds, mass, 'o-', ds, pred_mass, '*--'); +D = 15; +%plot(ds(1:D), mass(1:D), 'bo-', ds(1:D), pred_mass(1:D), 'g*--', ds(1:D), pred_mass2(1:D), 'k+-.'); +plot(ds(1:D), log(mass(1:D)), 'bo-', ds(1:D), log(pred_mass(1:D)), 'g*--', ds(1:D), log(pred_mass2(1:D)), 'k+-.'); + +grid on +xlabel('depth of hierarchy') +title('max num Q nodes in any clique vs. depth') +legend('actual', 'predicted') + +%previewfig(gcf, 'width', 3, 'height', 1.5, 'color', 'bw'); +%exportfig(gcf, '/home/cs/murphyk/WP/ConferencePapers/HHMM/clqsize2.eps', ... +% 'width', 3, 'height', 1.5, 'color', 'bw'); + +end + + +if 0 +for d=ds + effnumclqs(d) = length(find(Qsize(:,d)>0)); +end +ds = 1:10; +Qs = 2:10; +maxC = size(Qsize, 1); +cost = []; +cost_bound = []; +for qi=1:length(Qs) + Q = Qs(qi); + for d=ds + cost(d,qi) = 0; + for c=1:maxC + if length(Qsize(c,d) > 0) % this clique contains Q nodes + cost(d,qi) = cost(d,qi) + Q^Qsize(c,d)*2^Fsize(c,d); + end + end + %cost_bound(d,qi) = effnumclqs(d) * 8 * Q^(max(Qsize(:,d))); + cost_bound(d,qi) = (effnumclqs(d)*8) + Q^(max(Qsize(:,d))); + end +end + +qi=2; plot(ds, cost(:,qi), 'o-', ds, cost_bound(:,qi), '*--'); +end + + +if 0 +% convert numbers in cliques into names +for d=1:D + Fdecode(Fnodes(d)) = d; +end +for c=8:15 + clqs = J.cliques{c}; + fprintf('clique %d: ', c); + for k=clqs + if myismember(k, Qnodes) + fprintf('Q%d ', k) + elseif myismember(k, Fnodes) + fprintf('F%d ', Fdecode(k)) + elseif isequal(k, Onode) + fprintf('O ') + elseif myismember(k, Qnodes2) + fprintf('Q%d* ', k-ss) + else + error(['unrecognized node ' k]) + end + end + fprintf('\n'); +end +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm.m new file mode 100644 index 00000000..85ff7f6a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm.m @@ -0,0 +1,258 @@ +function [bnet, Qnodes, Fnodes, Onode] = mk_hhmm(varargin) +% MK_HHMM Make a Hierarchical HMM +% function [bnet, Qnodes, Fnodes, Onode] = mk_hhmm(...) +% +% e.g. 3-layer hierarchical HMM where level 1 only connects to level 2 +% and the parents of the observed node are levels 2 and 3. +% (This DBN is the same as Fig 10 in my tech report.) +% +% Q1 ----------> Q1 +% | \ ^ | +% | v / | +% | F2 ------/ | +% | ^ ^ \ | +% | / | \ | +% | / | || +% v | vv +% Q2----| --------> Q2 +% /| \ | ^| +% / | v | / | +% | | F3 --------/ | +% | | ^ \ | +% | v / v v +% | Q3 -----------> Q3 +% | | +% \ | +% v v +% O +% +% +% Optional arguments in name/value format [default value in brackets] +% +% Qsizes - sizes at each level [ none ] +% allQ - 1 means level i connects to all Q levels below, 0 means just to i+1 [0] +% transprob - transprob{d}(i,k,j) = P(Q(d,t)=j|Q(d,t-1)=i,Q(1:d-1,t)=k) ['leftright'] +% startprob - startprob{d}(k,j) = P(Q(d,t)=j|Q(1:d-1,t)=k) ['leftstart'] +% termprob - termprob{d}(k,j) = P(F(d,t)=2|Q(1:d-1,t)=k,Q(d,t)=j) for d>1 ['rightstop'] +% selfprop - prob of a self transition (termprob default = 1-selfprop) [0.8] +% Osize - size of O node +% discrete_obs - 1 means O is tabular_CPD, 0 means gaussian_CPD [0] +% Oargs - cell array of args to pass to the O CPD [ {} ] +% Ops - Q parents of O [Qnodes(end)] +% F1 - 1 means level 1 can finish (restart), else there is no F1->Q1 arc [0] +% clamp1 - 1 means we clamp the params of the Q nodes in slice 1 (Qt1params) [1] +% Note: the Qt1params are startprob, which should be shared with other slices. +% However, in the current implementation, the Qt1params will only be estimated +% from the initial state of each sequence. +% +% For d=1, startprob{1}(1,j) is only used in the first slice and +% termprob{1} is ignored, since we assume the top level never resets. +% Also, transprob{1}(i,j) can be used instead of transprob{1}(i,1,j). +% +% leftstart means the model always starts in state 1. +% rightstop means the model can only finish in its last state (Qsize(d)). +% unif means each state is equally like to reach any other +% rnd means the transition/starting probs are random (drawn from rand) +% +% Q1:QD in slice 1 are of type tabular_CPD +% Q1:QD in slice 2 are of type hhmmQ_CPD. +% F(2:D-1) is of type hhmmF_CPD, FD is of type tabular_CPD. + +args = varargin; +nargs = length(args); + +% get sizes of nodes and topology +Qsizes = []; +Osize = []; +allQ = 0; +Ops = []; +F1 = 0; +for i=1:2:nargs + switch args{i}, + case 'Qsizes', Qsizes = args{i+1}; + case 'Osize', Osize = args{i+1}; + case 'allQ', allQ = args{i+1}; + case 'Ops', Ops = args{i+1}; + case 'F1', F1 = args{i+1}; + end +end +if isempty(Qsizes), error('must specify Qsizes'); end +if Osize==0, error('must specify Osize'); end +D = length(Qsizes); +Qnodes = 1:D; + +if isempty(Ops), Ops = Qnodes(end); end + + +[intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, allQ, Ops, F1); +ss = length(intra); +names = {}; + +if F1 + Fnodes_ndx = Fnodes; +else + Fnodes_ndx = [-1 Fnodes]; % Fnodes(1) is a dummy index +end + +% set default params +discrete_obs = 0; +Oargs = {}; +startprob = cell(1,D); +startprob{1} = 'unif'; +for d=2:D + startprob{d} = 'leftstart'; +end +transprob = cell(1,D); +transprob{1} = 'unif'; +for d=2:D + transprob{d} = 'leftright'; +end +termprob = cell(1,D); +for d=2:D + termprob{d} = 'rightstop'; +end +selfprob = 0.8; +clamp1 = 1; + +for i=1:2:nargs + switch args{i}, + case 'discrete_obs', discrete_obs = args{i+1}; + case 'Oargs', Oargs = args{i+1}; + case 'startprob', startprob = args{i+1}; + case 'transprob', transprob = args{i+1}; + case 'termprob', termprob = args{i+1}; + case 'selfprob', selfprob = args{i+1}; + case 'clamp1', clamp1 = args{i+1}; + end +end + +ns = zeros(1,ss); +ns(Qnodes) = Qsizes; +ns(Onode) = Osize; +ns(Fnodes) = 2; + +dnodes = [Qnodes Fnodes]; +if discrete_obs + dnodes = [dnodes Onode]; +end +onodes = [Onode]; + +bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names); +eclass = bnet.equiv_class; + +for d=1:D + if d==1 + Qps = []; + elseif allQ + Qps = Qnodes(1:d-1); + else + Qps = Qnodes(d-1); + end + Qpsz = prod(ns(Qps)); + Qsz = ns(Qnodes(d)); + if isstr(startprob{d}) + switch startprob{d} + case 'unif', startprob{d} = mk_stochastic(ones(Qpsz, Qsz)); + case 'rnd', startprob{d} = mk_stochastic(rand(Qpsz, Qsz)); + case 'leftstart', startprob{d} = zeros(Qpsz, Qsz); startprob{d}(:,1) = 1; + end + end + if isstr(transprob{d}) + switch transprob{d} + case 'unif', transprob{d} = mk_stochastic(ones(Qsz, Qpsz, Qsz)); + case 'rnd', transprob{d} = mk_stochastic(rand(Qsz, Qpsz, Qsz)); + case 'leftright', + LR = mk_leftright_transmat(Qsz, selfprob); + temp = repmat(reshape(LR, [1 Qsz Qsz]), [Qpsz 1 1]); % transprob(k,i,j) + transprob{d} = permute(temp, [2 1 3]); % now transprob(i,k,j) + end + end + if isstr(termprob{d}) + switch termprob{d} + case 'unif', termprob{d} = mk_stochastic(ones(Qpsz, Qsz, 2)); + case 'rnd', termprob{d} = mk_stochastic(rand(Qpsz, Qsz, 2)); + case 'rightstop', + %termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i<Qsz + stopprob = 1-selfprob; + termprob{d} = zeros(Qpsz, Qsz, 2); + termprob{d}(:,Qsz,2) = stopprob; + termprob{d}(:,Qsz,1) = 1-stopprob; + termprob{d}(:,1:(Qsz-1),1) = 1; + otherwise, error(['unrecognized termprob ' termprob{d}]) + end + elseif d>1 % passed in termprob{d}(k,j) + temp = termprob{d}; + termprob{d} = zeros(Qpsz, Qsz, 2); + termprob{d}(:,:,2) = temp; + termprob{d}(:,:,1) = ones(Qpsz,Qsz) - temp; + end +end + + +% SLICE 1 + +for d=1:D + bnet.CPD{eclass(Qnodes(d),1)} = tabular_CPD(bnet, Qnodes(d), 'CPT', startprob{d}, 'adjustable', clamp1); +end + +if F1 + d = 1; + bnet.CPD{eclass(Fnodes_ndx(d),1)} = hhmmF_CPD(bnet, Fnodes_ndx(d), Qnodes(d), Fnodes_ndx(d+1), ... + 'termprob', termprob{d}); +end +for d=2:D-1 + if allQ + Qps = Qnodes(1:d-1); + else + Qps = Qnodes(d-1); + end + bnet.CPD{eclass(Fnodes_ndx(d),1)} = hhmmF_CPD(bnet, Fnodes_ndx(d), Qnodes(d), Fnodes_ndx(d+1), ... + 'Qps', Qps, 'termprob', termprob{d}); +end +bnet.CPD{eclass(Fnodes_ndx(D),1)} = tabular_CPD(bnet, Fnodes_ndx(D), 'CPT', termprob{D}); + +if discrete_obs + bnet.CPD{eclass(Onode,1)} = tabular_CPD(bnet, Onode, Oargs{:}); +else + bnet.CPD{eclass(Onode,1)} = gaussian_CPD(bnet, Onode, Oargs{:}); +end + +% SLICE 2 + +%for d=1:D +% bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, Qnodes, d, D, ... +% 'startprob', startprob{d}, 'transprob', transprob{d}, ... +% 'allQ', allQ); +%end + +d = 1; +if F1 + bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ... + 'Fbelow', Fnodes_ndx(d+1), ... + 'startprob', startprob{d}, 'transprob', transprob{d}); +else + bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, ... + 'Fbelow', Fnodes_ndx(d+1), ... + 'startprob', startprob{d}, 'transprob', transprob{d}); +end +for d=2:D-1 + if allQ + Qps = Qnodes(1:d-1); + else + Qps = Qnodes(d-1); + end + Qps = Qps + ss; % since all in slice 2 + bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ... + 'Fbelow', Fnodes_ndx(d+1), 'Qps', Qps, ... + 'startprob', startprob{d}, 'transprob', transprob{d}); +end +d = D; +if allQ + Qps = Qnodes(1:d-1); +else + Qps = Qnodes(d-1); +end +Qps = Qps + ss; % since all in slice 2 +bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ... + 'Qps', Qps, ... + 'startprob', startprob{d}, 'transprob', transprob{d}); diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo.m new file mode 100644 index 00000000..7a5fe57c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo.m @@ -0,0 +1,76 @@ +function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1) +% MK_HHMM_TOPO Make Hierarchical HMM topology +% function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1) +% +% D is the depth of the hierarchy +% If all_Q_to_Qs = 1, level i connects to all levels below, else just to i+1 [0] +% Ops are the Q parents of the observed node [Qnodes(end)] +% If F1=1, level 1 can finish (restart), else there is no F1->Q1 arc [0] + +Qnodes = 1:D; + +if nargin < 2, all_Q_to_Qs = 1; end +if nargin < 3, Ops = Qnodes(D); end +if nargin < 4, F1 = 0; end + +if F1 + Fnodes = 2*D:-1:D+1; % must number from bottom to top + Onode = 2*D+1; + ss = 2*D+1; +else + Fnodes = [-1 (2*D)-1:-1:D+1]; % Fnodes(1) is a dummy index + Onode = 2*D; + ss = 2*D; +end + +intra = zeros(ss); +intra(Ops, Onode) = 1; +for d=1:D-1 + if all_Q_to_Qs + intra(Qnodes(d), Qnodes(d+1:end)) = 1; + else + intra(Qnodes(d), Qnodes(d+1)) = 1; + end +end +for d=D:-1:3 + intra(Fnodes(d), Fnodes(d-1)) = 1; +end +if F1 + intra(Fnodes(2), Fnodes(1)) = 1; +end +if all_Q_to_Qs + if F1 + intra(Qnodes(1), Fnodes(1:end)) = 1; + else + intra(Qnodes(1), Fnodes(2:end)) = 1; + end + for d=2:D + intra(Qnodes(d), Fnodes(d:end)) = 1; + end +else + if F1 + intra(Qnodes(1), Fnodes([1 2])) = 1; + else + intra(Qnodes(1), Fnodes(2)) = 1; + end + for d=2:D-1 + intra(Qnodes(d), Fnodes([d d+1])) = 1; + end + intra(Qnodes(D), Fnodes(D)) = 1; +end + + +inter = zeros(ss); +for d=1:D + inter(Qnodes(d), Qnodes(d)) = 1; +end +if F1 + inter(Fnodes(1), Qnodes(1)) = 1; +end +for d=2:D + inter(Fnodes(d), Qnodes([d-1 d])) = 1; +end + +if ~F1 + Fnodes = Fnodes(2:end); % strip off dummy -1 term +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo_F1.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo_F1.m new file mode 100644 index 00000000..2fc5f912 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/mk_hhmm_topo_F1.m @@ -0,0 +1,65 @@ +function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo_F1(D, all_Q_to_Qs, Ops) +% MK_HHMM_TOPO Make Hierarchical HMM topology assuming level 1 can finish +% function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1) +% +% D is the depth of the hierarchy +% If all_Q_to_Qs = 1, level i connects to all levels below, else just to i+1 [0] +% Ops are the Q parents of the observed node [Qnodes(end)] +% If F1=1, level 1 can finish (restart), else there is no F1->Q1 arc [0] + +Qnodes = 1:D; + +if nargin < 2, all_Q_to_Qs = 1; end +if nargin < 3, Ops = Qnodes(D); end +if nargin < 4, F1 = 0; end + +if F1 + Fnodes = 2*D:-1:D+1; % must number from bottom to top + Onode = 2*D+1; + ss = 2*D+1; +else + Fnodes = (2*D)-1:-1:D+1; + Onode = 2*D; + ss = 2*D; +end + +intra = zeros(ss); +intra(Ops, Onode) = 1; +for d=1:D-1 + if all_Q_to_Qs + intra(Qnodes(d), Qnodes(d+1:end)) = 1; + else + intra(Qnodes(d), Qnodes(d+1)) = 1; + end +end +for d=D:-1:3 + intra(Fnodes(d), Fnodes(d-1)) = 1; +end +if F1 + intra(Fnodes(2), Fnodes(1)) = 1; +end +if all_Q_to_Qs + for d=1:D + intra(Qnodes(d), Fnodes(d:end)) = 1; + end +else + for d=1:D + if d < D + intra(Qnodes(d), Fnodes([d d+1])) = 1; + else + intra(Qnodes(d), Fnodes(d)) = 1; + end + end +end + +inter = zeros(ss); +for d=1:D + inter(Qnodes(d), Qnodes(d)) = 1; +end +for d=1:D + if d==1 + inter(Fnodes(d), Qnodes(d)) = 1; + else + inter(Fnodes(d), Qnodes([d-1 d])) = 1; + end +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/pretty_print_hhmm_parse.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/pretty_print_hhmm_parse.m new file mode 100644 index 00000000..81b965d5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/pretty_print_hhmm_parse.m @@ -0,0 +1,67 @@ +function pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, alphabet) +% function pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, alphabet) +% +% mpe(i,t) is the most probable value of node i at time t +% Qnodes(1:D), Fnodes = [F2 .. FD], Onode contain the node ids +% alphabet(i) is the i'th output symbol, or [] if don't want displayed + +T = size(mpe,2); +ncols = 20; +t1 = 1; t2 = min(T, t1+ncols-1); +while (t1 < T) + %fprintf('%d:%d\n', t1, t2); + if iscell(mpe) + print_block_cell(mpe(:,t1:t2), Qnodes, Fnodes, Onode, alphabet, t1); + else + print_block(mpe(:,t1:t2), Qnodes, Fnodes, Onode, alphabet, t1); + end + fprintf('\n\n'); + t1 = t2+1; t2 = min(T, t1+ncols-1); +end + +%%%%%% + +function print_block_cell(mpe, Qnodes, Fnodes, Onode, alphabet, start) + +D = length(Qnodes); +T = size(mpe, 2); +fprintf('%3d ', start:start+T-1); fprintf('\n'); +for d=1:D + for t=1:T + if (d > 1) & (mpe{Fnodes(d-1),t} == 2) + fprintf('%3d|', mpe{Qnodes(d), t}); + else + fprintf('%3d ', mpe{Qnodes(d), t}); + end + end + fprintf('\n'); +end +if ~isempty(alphabet) + a = cell2num(mpe(Onode,:)); + %fprintf('%3c ', alphabet(mpe{Onode,:})); + fprintf('%3c ', alphabet(a)) + fprintf('\n'); +end + + +%%%%%% + +function print_block(mpe, Qnodes, Fnodes, Onode, alphabet, start) + +D = length(Qnodes); +T = size(mpe, 2); +fprintf('%3d ', start:start+T-1); fprintf('\n'); +for d=1:D + for t=1:T + if (d > 1) & (mpe(Fnodes(d-1),t) == 2) + fprintf('%3d|', mpe(Qnodes(d), t)); + else + fprintf('%3d ', mpe(Qnodes(d), t)); + end + end + fprintf('\n'); +end +if ~isempty(alphabet) + fprintf('%3c ', alphabet(mpe(Onode,:))); + fprintf('\n'); +end diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/remove_hhmm_end_state.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/remove_hhmm_end_state.m new file mode 100644 index 00000000..1ef9ded9 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/remove_hhmm_end_state.m @@ -0,0 +1,41 @@ +function [transprob, termprob] = remove_hhmm_end_state(A) +% REMOVE_END_STATE Infer transition and termination probabilities from automaton with an end state +% [transprob, termprob] = remove_end_state(A) +% +% A(i,k,j) = Pr( i->j | Qps=k), where i in 1:Q, j in 1:(Q+1), and Q+1 is the end state +% This implements the equation in footnote 3 of my NIPS 01 paper, +% transprob(i,k,j) = \tilde{A}_k(i,j) +% termprob(k,j) = \tau_k(j) +% +% For the top level, the k index is missing. + +Q = size(A,1); +toplevel = (ndims(A)==2); +if toplevel + Qk = 1; + A = reshape(A, [Q 1 Q+1]); +else + Qk = size(A, 2); +end + +transprob = A(:, :, 1:Q); +term = A(:,:,Q+1)'; % term(k,j) = P(Qj -> end | k) +termprob = term; +%termprob = zeros(Qk, Q, 2); +%termprob(:,:,2) = term; +%termprob(:,:,1) = 1-term; + +for k=1:Qk + for i=1:Q + for j=1:Q + denom = (1-termprob(k,i)); + denom = denom + (denom==0)*eps; + transprob(i,k,j) = transprob(i,k,j) / denom; + end + end +end + +if toplevel + termprob = squeeze(termprob); + transprob = squeeze(transprob); +end |