Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/general/mk_dbn.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function bnet = mk_dbn(intra, inter, node_sizes, varargin)
wolffd@0 2 % MK_DBN Make a Dynamic Bayesian Network.
wolffd@0 3 %
wolffd@0 4 % BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs
wolffd@0 5 % from i in slice t to j in slice t iff intra(i,j) = 1, and
wolffd@0 6 % from i in slice t to j in slice t+1 iff inter(i,j) = 1,
wolffd@0 7 % for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1.
wolffd@0 8 % node_sizes(i) is the number of values node i can take on.
wolffd@0 9 % The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
wolffd@0 10 % See also mk_bnet.
wolffd@0 11 %
wolffd@0 12 % Optional arguments [default in brackets]
wolffd@0 13 % 'discrete' - list of discrete nodes [1:n]
wolffd@0 14 % 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ]
wolffd@0 15 % 'eclass1' - equiv class for slice 1 [1:n]
wolffd@0 16 % 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]
wolffd@0 17 % equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j},
wolffd@0 18 % i.e., nodes i and j have tied parameters.
wolffd@0 19 % 'intra1' - topology of first slice, if different from others
wolffd@0 20 % 'names' - a cell array of strings to be associated with nodes 1:n [{}]
wolffd@0 21 % This creates an associative array, so you write e.g.
wolffd@0 22 % 'evidence(bnet.names{'bar'}) = 42' instead of 'evidence(2} = 42'
wolffd@0 23 % assuming names = { 'foo', 'bar', ...}.
wolffd@0 24 %
wolffd@0 25 % For backwards compatibility with BNT2, arguments can also be specified as follows
wolffd@0 26 % bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1)
wolffd@0 27 %
wolffd@0 28 % After calling this function, you must specify the parameters (conditional probability
wolffd@0 29 % distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc.
wolffd@0 30
wolffd@0 31
wolffd@0 32 n = length(intra);
wolffd@0 33 ss = n;
wolffd@0 34 bnet.nnodes_per_slice = ss;
wolffd@0 35 bnet.intra = intra;
wolffd@0 36 bnet.inter = inter;
wolffd@0 37 bnet.intra1 = intra;
wolffd@0 38 dag = zeros(2*n);
wolffd@0 39 dag(1:n,1:n) = bnet.intra1;
wolffd@0 40 dag(1:n,(1:n)+n) = bnet.inter;
wolffd@0 41 dag((1:n)+n,(1:n)+n) = bnet.intra;
wolffd@0 42 bnet.dag = dag;
wolffd@0 43 bnet.names = {};
wolffd@0 44
wolffd@0 45 directed = 1;
wolffd@0 46 if ~acyclic(dag,directed)
wolffd@0 47 error('graph must be acyclic')
wolffd@0 48 end
wolffd@0 49
wolffd@0 50
wolffd@0 51 bnet.eclass1 = 1:n;
wolffd@0 52 %bnet.eclass2 = (1:n)+n;
wolffd@0 53 bnet.eclass2 = bnet.eclass1;
wolffd@0 54 for i=1:ss
wolffd@0 55 if isequal(parents(dag, i+ss), parents(dag, i)+ss)
wolffd@0 56 %fprintf('%d has isomorphic parents, eclass %d\n', i, bnet.eclass2(i))
wolffd@0 57 else
wolffd@0 58 bnet.eclass2(i) = max(bnet.eclass2) + 1;
wolffd@0 59 %fprintf('%d has non isomorphic parents, eclass %d\n', i, bnet.eclass2(i))
wolffd@0 60 end
wolffd@0 61 end
wolffd@0 62
wolffd@0 63 dnodes = 1:n;
wolffd@0 64 bnet.observed = [];
wolffd@0 65
wolffd@0 66 if nargin >= 4
wolffd@0 67 args = varargin;
wolffd@0 68 nargs = length(args);
wolffd@0 69 if ~isstr(args{1})
wolffd@0 70 if nargs >= 1, dnodes = args{1}; end
wolffd@0 71 if nargs >= 2, bnet.eclass1 = args{2}; end
wolffd@0 72 if nargs >= 3, bnet.eclass2 = args{3}; end
wolffd@0 73 if nargs >= 4, bnet.intra1 = args{4}; end
wolffd@0 74 else
wolffd@0 75 for i=1:2:nargs
wolffd@0 76 switch args{i},
wolffd@0 77 case 'discrete', dnodes = args{i+1};
wolffd@0 78 case 'observed', bnet.observed = args{i+1};
wolffd@0 79 case 'eclass1', bnet.eclass1 = args{i+1};
wolffd@0 80 case 'eclass2', bnet.eclass2 = args{i+1};
wolffd@0 81 case 'intra1', bnet.intra1 = args{i+1};
wolffd@0 82 %case 'ar_hmm', bnet.ar_hmm = args{i+1}; % should check topology
wolffd@0 83 case 'names', bnet.names = assocarray(args{i+1}, num2cell(1:n));
wolffd@0 84 otherwise,
wolffd@0 85 error(['invalid argument name ' args{i}]);
wolffd@0 86 end
wolffd@0 87 end
wolffd@0 88 end
wolffd@0 89 end
wolffd@0 90
wolffd@0 91
wolffd@0 92 bnet.observed = sort(bnet.observed); % for comparing sets
wolffd@0 93 ns = node_sizes;
wolffd@0 94 bnet.node_sizes_slice = ns(:)';
wolffd@0 95 bnet.node_sizes = [ns(:) ns(:)];
wolffd@0 96
wolffd@0 97 cnodes = mysetdiff(1:n, dnodes);
wolffd@0 98 bnet.dnodes_slice = dnodes;
wolffd@0 99 bnet.cnodes_slice = cnodes;
wolffd@0 100 bnet.dnodes = [dnodes dnodes+n];
wolffd@0 101 bnet.cnodes = [cnodes cnodes+n];
wolffd@0 102
wolffd@0 103 bnet.equiv_class = [bnet.eclass1(:) bnet.eclass2(:)];
wolffd@0 104 bnet.CPD = cell(1,max(bnet.equiv_class(:)));
wolffd@0 105 eclass = bnet.equiv_class(:);
wolffd@0 106 E = max(eclass);
wolffd@0 107 bnet.rep_of_eclass = zeros(1,E);
wolffd@0 108 for e=1:E
wolffd@0 109 mems = find(eclass==e);
wolffd@0 110 bnet.rep_of_eclass(e) = mems(1);
wolffd@0 111 end
wolffd@0 112
wolffd@0 113 ss = n;
wolffd@0 114 onodes = bnet.observed;
wolffd@0 115 hnodes = mysetdiff(1:ss, onodes);
wolffd@0 116 bnet.hidden_bitv = zeros(1,2*ss);
wolffd@0 117 bnet.hidden_bitv(hnodes) = 1;
wolffd@0 118 bnet.hidden_bitv(hnodes+ss) = 1;
wolffd@0 119
wolffd@0 120 bnet.parents = cell(1, 2*ss);
wolffd@0 121 for i=1:ss
wolffd@0 122 bnet.parents{i} = parents(bnet.dag, i);
wolffd@0 123 bnet.parents{i+ss} = parents(bnet.dag, i+ss);
wolffd@0 124 end
wolffd@0 125
wolffd@0 126 bnet.auto_regressive = zeros(1,ss);
wolffd@0 127 % ar(i)=1 means (observed) node i depends on i in the previous slice
wolffd@0 128 for o=bnet.observed(:)'
wolffd@0 129 if any(bnet.parents{o+ss} <= ss)
wolffd@0 130 bnet.auto_regressive(o) = 1;
wolffd@0 131 end
wolffd@0 132 end
wolffd@0 133