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