function pot = marginal_nodes(bnet, nodes, t, c)
% MARGINAL_NODES Compute the marginal distribution Pr(X|evidence) on the specified set of nodes X.
%
% M = MARGINAL_NODES(BNET, NODES) computes the marginal distribution on 'nodes'.
% M = MARGINAL_NODES(BNET, NODES, T) compute the marginal distribution on 'nodes' in slice t
% (and possibly t+1, if any(nodes > ss), where ss = slice size).
% M = MARGINAL_NODES(BNET, NODES, T, C) specifies that clique C contains 'nodes'.
% (This is used by MARGINAL_FAMILY if algo=='jtree'.)
%
% 'm' is a structure; m.T is a table with the distribution over the discrete nodes (if any).
% m.mu and m.Sigma contains the first and second moments of the distribution over the cts nodes
% (if any).
%
% In general, it is not possible to compute the posterior on an arbitrary subset of nodes;
% it depends on which inference algorithm you are using. If bnet.algo=='jtree' or 'bk', you can
% call SCHEDULE_MARGINALS before you call ENTER_EVIDENCE to ensure that there will be a clique
% containing the marginal you require. If bnet=='hmm' or 'lds', you can get the marginal on any
% subset of nodes in adjacent time-slices. For sampling and Pearl, you can only get marginals
% on subsets of families (see MARGINAL_FAMILY).

%if ~exist('t'), t = 1; end
%if ~exist('c'), c = -1; end
if nargin < 3, t = 1; end
if nargin < 4, c = -1; end

ss = bnet.slice_size;

switch bnet.algo
  case {'hmm', 'lds'}
    if all(nodes<=ss) % all fit within one slice
      bigm = bnet.hmm.one_slice_marginal{t};
    else
      bigm = bnet.hmm.two_slice_marginal{t};
    end
  case 'bk',
    % post{1} contains slices 1 and 2, post{2} 2 and 3, ..., post{T-1} T-1 and T
    % If we are only doing filtering, and we want the marginal on a node in slice t,
    % we should avoid including the evidence from slice t+1, so we use post{t-1} instead of post{t}
    % For slice 1, we should use post1 instead of post{1}.
    if t==1 & subset(nodes, (1:ss))
      if bnet.filter_only
	c = find_clique(bnet.jtree1, nodes);
	bigm = bnet.dbn.post1{c};
      else
	c = find_clique(bnet.jtree, nodes);
	bigm = bnet.dbn.post{1}{c};
      end
    else
      if subset(nodes, 1:ss)
	c = find_clique(bnet.jtree, nodes+ss); % find clique in slice 2
	bigm = bnet.dbn.post{t-1}{c};
      else
	c = find_clique(bnet.jtree, nodes);
	bigm = bnet.dbn.post{t}{c};
      end
    end
  case 'jtree',
    if c == -1
      if length(nodes)==1
	c = find_clique_containing_family(bnet, nodes(1), t);
      else
	c = find_clique(bnet.jtree, nodes + (t-1)*ss);
      end
    end
    bigm = bnet.jtree.clpot{c};
  case {'ls', 'lw'}
    bigm = marginal_from_samples(bnet.sample.samples, bnet.sample.weights, nodes + (t-1)*ss, ...
	bnet.node_sizes, bnet.discrete_nodes, bnet.cts_nodes);
  case 'sof'
    if any(nodes>ss) % spans slices t and t+1
      nsamples = min(size(bnet.sof.samples{t}, 1), size(bnet.sof.samples{t+1}, 1));
      samples = blockmatrix(nsamples, [], 1, [bnet.blocks_slice bnet.blocks_slice]);
      %%%samples = blockmatrix(nsamples, [bnet.blocks_slice bnet.blocks_slice]);
      samples(:,1:ss) = bnet.sof.samples{t}(1:nsamples,:);
      samples(:,(1:ss)+ss) = bnet.sof.samples{t+1}(1:nsamples, :);
      weights = bnet.sof.weights{t}(1:nsamples); % what about the other weights?
    else
      nsamples = size(bnet.sof.samples{t}, 1);
      samples = blockmatrix(nsamples, [], 1, bnet.blocks_slice);
      %%%samples = blockmatrix(nsamples, bnet.blocks_slice);
      samples = bnet.sof.samples{t};
      weights = bnet.sof.weights{t};
    end
    bigm = marginal_from_samples(samples, weights, nodes, bnet.node_sizes, bnet.discrete_nodes,  bnet.cts_nodes);
    bigm = set_timeslice_belief(bigm, t, bnet);
  case 'pearl',
    assert(t==1);
    if length(nodes)==1
      bigm = bnet.pearl.marginal{nodes(1)};
    else
      error('can''t use marginal_nodes with pearl');
    end
  case 'cutset',
    assert(t==1);
    if length(nodes)==1
      bigm = bnet.cutset.marginal{i};
    else
      error('can''t use marginal_nodes with cutset');
    end
  otherwise,
    error(['can''t use ' bnet.algo ' with marginal_nodes']);
end

% change to global coords before marginalizing bigm
nodes = nodes + (t-1)*ss;
bigm = ensure_moment_form(bigm);
m = marginalize_potential(bigm, nodes);


% Convert CG potentials to simpler form if possible.
switch m.type
  case 'cg',
    if subset(nodes, bnet.discrete_nodes)
      m = cg_to_table(m);
    elseif subset(nodes, bnet.cts_nodes)
      m = cg_to_gpot(m);
    end
end

pot = add_evidence_to_marginal(bnet, m);



