matthiasm@8: function [marginal, loglik] = marginal_nodes(engine, query)
matthiasm@8: % MARGINAL_NODES Compute the marginal on the specified query nodes (enumerative_inf)
matthiasm@8: % [marginal, loglik] = marginal_nodes(engine, query)
matthiasm@8: 
matthiasm@8: 
matthiasm@8: if isempty(query) & nargout < 2
matthiasm@8:   marginal.T = 1;
matthiasm@8:   marginal.domain = [];
matthiasm@8:   return;
matthiasm@8: end
matthiasm@8: 
matthiasm@8: evidence = engine.evidence;
matthiasm@8: bnet = bnet_from_engine(engine);
matthiasm@8: assert(isempty(bnet.cnodes));
matthiasm@8: n = length(bnet.dag);
matthiasm@8: observed = ~isemptycell(evidence);
matthiasm@8: vals = cat(1,evidence{observed});
matthiasm@8: vals = vals(:)';
matthiasm@8: ns = bnet.node_sizes;
matthiasm@8: 
matthiasm@8: sz = ns(query);
matthiasm@8: T = 0*myones(sz);
matthiasm@8: p = 0;
matthiasm@8: for i=1:prod(ns)
matthiasm@8:   inst = ind2subv(ns, i); % i'th instantiation
matthiasm@8:   if isempty(vals) | inst(observed) == vals % agrees with evidence
matthiasm@8:     prob = exp(log_lik_complete(bnet, num2cell(inst(:))));
matthiasm@8:     p = p + prob;
matthiasm@8:     v = inst(query);
matthiasm@8:     j = subv2ind(sz, v);
matthiasm@8:     T(j) = T(j) + prob;
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: [T, lik] = normalise(T);
matthiasm@8: lik = p;
matthiasm@8: loglik = log(lik);
matthiasm@8: 
matthiasm@8: Tsmall = shrink_obs_dims_in_table(T, query, evidence);
matthiasm@8: marginal.domain = query;
matthiasm@8: marginal.T = Tsmall;