Daniel@0: function marginal = marginal_nodes(engine, nodes, t, fam)
Daniel@0: % MARGINAL_NODES Compute the marginal on the specified query nodes (bk)
Daniel@0: %
Daniel@0: %   marginal = marginal_nodes(engine, i, t)
Daniel@0: % returns Pr(X(i,t) | Y(1:T)), where X(i,t) is the i'th node in the t'th slice.
Daniel@0: % If enter_evidence used filtering instead of smoothing, this will return  Pr(X(i,t) | Y(1:t)).
Daniel@0: %
Daniel@0: %   marginal = marginal_nodes(engine, query, t)
Daniel@0: % returns Pr(X(query(1),t), ... X(query(end),t) | Y(1:T)),
Daniel@0: % where X(q,t) is the q'th node in the t'th slice. If q > ss (slice size), this is equal
Daniel@0: % to X(q mod ss, t+1). That is, 't' specifies the time slice of the earliest node.
Daniel@0: % 'query' cannot span more than 2 time slices.
Daniel@0: % Example:
Daniel@0: % Consider a DBN with 2 nodes per slice.
Daniel@0: % Then t=2, nodes=[1 3] refers to node 1 in slice 2 and node 1 in slice 3.
Daniel@0: 
Daniel@0: if nargin < 3, t = 1; end
Daniel@0: if nargin < 4, fam = 0; else fam = 1; end
Daniel@0: 
Daniel@0: 
Daniel@0: % clpot{t} contains slice t-1 and t
Daniel@0: % Example
Daniel@0: % clpot #: 1    2    3
Daniel@0: % slices:  1  1,2  2,3
Daniel@0: % For filtering, we must take care not to take future evidence into account.
Daniel@0: % For smoothing, clpot{1} does not exist.
Daniel@0: 
Daniel@0: bnet = bnet_from_engine(engine);
Daniel@0: ss = length(bnet.intra);
Daniel@0: 
Daniel@0: nodes2 = nodes;
Daniel@0: if ~engine.filter
Daniel@0:   if t < engine.T
Daniel@0:     slice = t+1;
Daniel@0:   else % earliest t is T, so all nodes fit in one slice
Daniel@0:     slice = engine.T;
Daniel@0:     nodes2 = nodes + ss;
Daniel@0:   end
Daniel@0: else
Daniel@0:   if t == 1
Daniel@0:    slice = 1;
Daniel@0:   else
Daniel@0:     if all(nodes<=ss)
Daniel@0:       slice = t;
Daniel@0:       nodes2 = nodes + ss;
Daniel@0:     elseif t == engine.T
Daniel@0:       slice = t;
Daniel@0:     else
Daniel@0:       slice = t + 1;
Daniel@0:     end
Daniel@0:   end
Daniel@0: end
Daniel@0:   
Daniel@0: if engine.filter & t==1
Daniel@0:   c = clq_containing_nodes(engine.sub_engine1, nodes2, fam);
Daniel@0: else
Daniel@0:   c = clq_containing_nodes(engine.sub_engine, nodes2, fam);
Daniel@0: end
Daniel@0: assert(c >= 1);
Daniel@0: bigpot = engine.clpot{c, slice};
Daniel@0: 
Daniel@0: pot = marginalize_pot(bigpot, nodes2);
Daniel@0: marginal = pot_to_marginal(pot);
Daniel@0: 
Daniel@0: % we convert the domain to the unrolled numbering system
Daniel@0: % so that update_ess extracts the right evidence.
Daniel@0: marginal.domain = nodes+(t-1)*ss;