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