function [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, dom, sizes)
% MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian.
% [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, dom, sizes)
% 
% sizes(i) is the size of the i'th block in domain.
% logp = log (weight * normal_coef(Sigma))
% logpX = log (weight * normal_coef(SigmaX))
% where normal_coef(S) is the normalizing constant for a Gaussian with cov. S
% and weight is an arbitrary positive constant.

if isempty(dom)
  assert(isempty(X));
  assert(isempty(Y));
else
  assert(isequal(dom, union(X, Y))); % can fail if all args are empty
end

if isempty(Y)
  logpX = logp;
  muX = mu;
  SXX = Sigma;
  return;
end

[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, dom, sizes);

%logpX = logp - log(normal_coef(Sigma)) + log(normal_coef(SXX)); % COEF
logpX = logp; % Lauritzen (1996) p161
