diff toolboxes/FullBNT-1.0.7/bnt/general/convert_dbn_CPDs_to_tables.m @ 0:e9a9cd732c1e tip

first hg version after svn

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolboxes/FullBNT-1.0.7/bnt/general/convert_dbn_CPDs_to_tables.m	Tue Feb 10 15:05:51 2015 +0000
@@ -0,0 +1,201 @@
+function CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
+% CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables
+% CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
+%
+% CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
+% All hidden nodes are assumed to be discrete.
+% We assume the observed nodes are the same in every slice.
+%
+% Evaluating the conditional likelihood of long evidence sequences can be very slow,
+% so we take pains to vectorize where possible.
+
+[ss T] = size(evidence);
+%obs_bitv = ~isemptycell(evidence(:));
+obs_bitv = zeros(1, 2*ss);
+obs_bitv(bnet.observed) = 1;
+obs_bitv(bnet.observed+ss) = 1;
+
+ns = bnet.node_sizes(:);
+CPDpot = cell(ss,T); 
+
+for n=1:ss
+  % slice 1
+  t = 1;
+  ps = parents(bnet.dag, n);
+  e = bnet.equiv_class(n, 1);
+  if ~any(obs_bitv(ps))
+    CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
+  else
+    CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
+  end
+
+% special cases: c=child, p=parents, d=discrete, h=hidden, 1sl=1slice
+% if c=h=1 then c=d=1, since hidden nodes must be discrete
+% c=h c=d p=h p=d 1sl method
+% ---------------------------
+% 1   1   1   1   -   replicate CPT
+% -   1   -   1   -   evaluate CPT on evidence *
+% 0   1   1   1   1   dhmm
+% 0   0   1   1   1   ghmm
+% other               loop
+%
+% * = any subset of the domain may be observed
+
+% Example where all of the special cases occur - a hierarchical HMM
+% where the top layer (G) and leaves (Y) are observed and
+% all nodes are discrete except Y.
+% (O turns on if Y is an outlier)
+
+% G ---------> G 
+% |            |
+% v            v
+% S  --------> S
+% |            |
+% v            v
+% Y            Y
+% ^            ^
+% |            |
+% O            O
+
+% Evaluating P(yt|St,Ot) is the ghmm case
+% Evaluating P(St|S(t-1),gt) is the eval CPT case
+% Evaluating P(gt|g(t-1) is the eval CPT case (hdom = [])
+% Evaluating P(Ot) is the replicated CPT case
+
+% Cts parents (e.g., inputs) would require an additional special case for speed
+
+
+  % slices 2..T
+  [ss T] = size(evidence);
+  self = n+ss;
+  ps = parents(bnet.dag, self);
+  e = bnet.equiv_class(n, 2);
+
+  if 1
+  debug = 0;
+  hidden_child = ~obs_bitv(n);
+  discrete_child = myismember(n, bnet.dnodes);
+  hidden_ps = all(~obs_bitv(ps));
+  discrete_ps = mysubset(ps, bnet.dnodes);
+  parents_in_same_slice = all(ps > ss);
+  
+  if hidden_child & discrete_child & hidden_ps & discrete_ps
+    CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  elseif discrete_child & discrete_ps
+    CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  elseif discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
+    CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  elseif ~discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
+    CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
+  else
+    if debug, fprintf('node %d, slow\n', n); end
+    for t=2:T
+      CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    end
+  end
+  end
+  
+  if 0
+  for t=2:T
+    CPDpot2{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
+    if ~approxeq(CPDpot{n,t}, CPDpot2{n,t})
+      fprintf('CPDpot n=%d, t=%d\n',n,t);
+      keyboard
+    end
+  end
+  end
+
+  
+end
+
+
+
+
+%%%%%%%
+function CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+[ss T] = size(evidence);
+if debug, fprintf('node %d, repl\n', n); end
+e = bnet.equiv_class(n, 2);
+CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []);
+CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]);
+
+
+
+%%%%%%%
+function CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+[ss T] = size(evidence);
+self = n+ss;
+ps = parents(bnet.dag, self);
+e = bnet.equiv_class(n, 2);
+ns = bnet.node_sizes(:);
+% Example: given CPT(p1, p2, p3, p4, c), where p1,p3 are observed
+% we create CPT([p2 p4 c], [p1 p3]).
+% We then convert all observed p1,p3 into indices ndx
+% and return CPT(:, ndx)
+CPT = CPD_to_CPT(bnet.CPD{e});
+domain = [ps self];
+% if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2,
+% odom = [3 8], hdom = 7
+odom_rel = find(obs_bitv(domain));
+hdom_rel = find(~obs_bitv(domain));
+odom = domain(odom_rel);
+hdom = domain(hdom_rel);
+if isempty(hdom)
+  CPT = CPT(:);
+else
+  CPT = permute(CPT, [hdom_rel odom_rel]);
+  CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom)));
+end
+parents_in_same_slice = all(ps > ss);
+if parents_in_same_slice
+  if debug, fprintf('node %d eval 1 slice\n', n); end
+  data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1
+else
+  if debug, fprintf('node %d eval 2 slice\n', n); end
+  % there's probably a way of vectorizing this...
+  data = zeros(length(odom), T-1);
+  for t=2:T
+    ev = evidence(:,t-1:t);
+    ev = ev(:);
+    ev2 = ev(odom);
+    data(:,t-1) = cat(1, ev2{:});
+    %data(:,t-1) = cell2num(ev2);
+  end
+end
+ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t)
+if isempty(hdom)
+  CPDpot(n,2:T) = num2cell(CPT(ndx)); % a cell array of floats
+else
+  CPDpot(n,2:T) = num2cell(CPT(:, ndx), 1); % a cell array of column vectors
+end
+
+%%%%%%%
+function CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+if debug, fprintf('node %d, dhmm\n', n); end
+[ss T] = size(evidence);
+self = n+ss;
+ps = parents(bnet.dag, self);
+e = bnet.equiv_class(n, 2);
+ns = bnet.node_sizes(:);
+CPT = CPD_to_CPT(bnet.CPD{e});
+CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents?
+%obslik = mk_dhmm_obs_lik(cell2num(evidence(n,2:T)), CPT);
+obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT);
+CPDpot(n,2:T) = num2cell(obslik, 1);
+
+
+%%%%%%%
+function CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug)
+
+if debug, fprintf('node %d, ghmm\n', n); end
+[ss T] = size(evidence);
+e = bnet.equiv_class(n, 2);
+S = struct(bnet.CPD{e}); 
+ev2 = cell2num(evidence(n,2:T));
+%obslik = mk_ghmm_obs_lik(ev2, S.mean, S.cov);
+obslik = eval_pdf_cond_gauss(ev2, S.mean, S.cov);
+CPDpot(n,2:T) = num2cell(obslik, 1);
+