enter_evidence_debug.m

来自「Bayes网络工具箱」· M 代码 · 共 167 行

function [engine, loglik] = enter_evidence(engine, evidence, filter)
% ENTER_EVIDENCE Add the specified evidence to the network (frontier)
% [engine, loglik] = enter_evidence(engine, evidence, filter)
%
% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
% If filter = 1, we do filtering, otherwise smoothing (default).

if nargin < 3, filter = 0; end

[ss T] = size(evidence);
bnet = bnet_from_engine(engine);
ns = repmat(bnet.node_sizes_slice(:), 1, T);
eclass = [bnet.equiv_class(:,1) repmat(bnet.equiv_class(:,2), 1, T-1)];
onodes = find(~isemptycell(evidence));
cnodes = unroll_set(bnet.cnodes(:), ss, T);
big_dag = unroll_dbn_topology(bnet.intra1, bnet.intra, bnet.inter, T);
pot_type = determine_pot_type(onodes, cnodes, big_dag);

ens = ns;
ens(onodes) = 1;
ns = ns(:)';
cnodes = cnodes(:)';

% Convert evidence-specific CPDs to potentials
CPD = cell(ss,T);
for t=1:T
  for i=1:ss
    fam = family(bnet.dag, i, t);
    CPD{i,t} = CPD_to_pot(pot_type, bnet.CPD{eclass(i,t)}, fam, ns, cnodes, evidence);   
  end
end


% FORWARDS
fwd = cell(ss,T);
ll = zeros(1,T);
S = 2*ss; % num. intermediate frontiers to get from t to t+1
frontier = cell(S,T);

% Start with empty frontier, and add each node in slice 1
init = mk_initial_pot(pot_type, [], ns, cnodes, onodes);  
t = 1;
s = 1;
j = 1;
[frontier{s,t}, fndx1{s,t}, CPDndx1{s,t}] = update(init, ens, j, 1, CPD{j}, engine.fdom1{s}, pot_type, ns, cnodes, onodes);
fwd{j} = frontier{s,t};
for s=2:ss
  j = s; % add node j at step s
  [frontier{s,t}, fndx1{s,t}, CPDndx1{s,t}] = update(frontier{s-1,t}, ens, j, 1, CPD{j}, engine.fdom1{s}, pot_type, ns, cnodes, onodes);
  fwd{j} = frontier{s,t};
end
frontier{S,t} = frontier{ss,t};
[frontier{S,t}, ll(1)] = normalize_pot(frontier{S,t});

% Now move frontier from slice to slice
OPS = engine.ops;
add = OPS>0;
nodes = [zeros(S,1) unroll_set(abs(OPS(:)), ss, T-1)];
for t=2:T
  offset = (t-2)*ss;
  for s=1:S
    if s==1
      prev_ndx = (t-2)*S + S; % S,t-1
    else
      prev_ndx = (t-1)*S + s-1; % s-1,t
    end
    j = nodes(s,t);
    [frontier{s,t}, fndx{s,t}, CPDndx{s,t}] = update(frontier{prev_ndx}, ens, j, add(s), CPD{j}, engine.fdom{s}+offset, pot_type, ns, cnodes, onodes);
    if add(s)
      fwd{j} = frontier{s,t};
    end
  end
  [frontier{S,t}, ll(t)] = normalize_pot(frontier{S,t});
end
loglik = sum(ll);

engine.fwd = fwd;
engine.fndx = fndx;
engine.CPDndx = CPDndx;
engine.fndx1 = fndx1;
engine.CPDndx1 = CPDndx1;

if 1 % filter
  engine.fwdback = fwd;
  return;
end


% BACKWARDS
back = cell(ss,T);
add = ~add; % forwards add = backwards remove 
frontier = cell(S,T+1);
t = T;
dom = (1:ss) + (t-1)*ss;
frontier{1,T+1} = mk_initial_pot(pot_type, dom, ns, cnodes, onodes); % all 1s for last slice
for t=T:-1:2
  offset = (t-2)*ss;
  for s=S:-1:1 % reverse order
    if s==S
      prev_ndx = t*S + 1; % 1,t+1
    else
      prev_ndx = (t-1)*S + (s+1); % s+1,t
    end
    j = nodes(s,t);
    if ~add(s)
      back{j} = frontier{prev_ndx}; % save frontier before removing
    end
    frontier{s,t} = rev_update(frontier{prev_ndx}, j, add(s), CPD{j}, engine.fdom{s}+offset, pot_type, ns, cnodes, onodes);
  end
  frontier{1,t} = normalize_pot(frontier{1,t});
end
% Remove each node in first slice until left with empty set
t = 1;
frontier{ss+1,t} = frontier{1,2};
for s=ss:-1:1
  j = s; % remove node j at step s
  back{j} = frontier{s+1,t};
  frontier{s,t} = rev_update(frontier{s+1,t}, j, 0, CPD{j}, 1:s, pot_type, ns, cnodes, onodes);
end

% COMBINE
for t=1:T
  for i=1:ss
   engine.fwdback{i,t} = normalize_pot(multiply_pots(fwd{i,t}, back{i,t}));
  end
end

% for debugging
engine.back = back;


%%%%%%%%%%
function [new_frontier, fndx, CPDndx]  = update(old_frontier, ens, j, add, CPD, newdom, pot_type, ns, cnodes, onodes)

olddom = domain_pot(old_frontier);
if add
  new_frontier = mk_initial_pot(pot_type, newdom, ns, cnodes, onodes);      
  new_frontier = multiply_by_pot(new_frontier, old_frontier);
  new_frontier = multiply_by_pot(new_frontier, CPD);

  fndx = mk_multiply_table_ndx(newdom, olddom, ens);
  temp = struct(CPD);
  CPDndx = mk_multiply_table_ndx(newdom, temp.domain, ens);
else
  new_frontier = marginalize_pot(old_frontier, mysetdiff(domain_pot(old_frontier), j));    
  
  fndx = mk_marginalise_table_ndx(olddom, newdom, ens);
  CPDndx = [];
end


%%%%%%
function new_frontier = rev_update(old_frontier, j, add, CPD, newdom, pot_type, ns, cnodes, onodes)

if add
  % add: extend domain to include i by multiplying by 1
  new_frontier = mk_initial_pot(pot_type, newdom, ns, cnodes, onodes);      
  new_frontier = multiply_by_pot(new_frontier, old_frontier);
else 
  % remove: multiply in CPT and then marginalize it out
  temp = mk_initial_pot(pot_type, newdom, ns, cnodes, onodes);      
  temp = multiply_by_pot(temp, old_frontier);
  temp = multiply_by_pot(temp, CPD);
  new_frontier = marginalize_pot(temp, mysetdiff(domain_pot(temp), j));
end