function [mu2, Sigma2, ll] = condition_gaussian(mu, Sigma, X, Y, dom, y, sizes)
% CONDITION_GAUSSIAN Compute Pr(X|Y=y) where X and Y are jointly Gaussian.
% [mu2, Sigma2, ll] = condition_gaussian(mu, Sigma, X, Y, dom, y, sizes)
%
% dom = union(X,Y) and sizes(i) is the size of the i'th block.
% mu2 = E[X|Y=y] is a n*1 vector and Sigma2 = Cov[X|Y=y] is a n*n matrix, where n = sum(sizes(X)).
% ll = log Pr(Y=y)
% e.g., conditional_gaussian(mu, Sigma, [1 3], [2], [1 2 3], [-0.5 0.5 2.1], [2 3 4]) so n=6.

if isempty(y)
  mu2 = mu;
  Sigma2 = Sigma;
  ll = 0;
  return;
end

y = y(:);

if length(Y)==length(dom) % instantiating every variable
  mu2 = y;
  Sigma2 = zeros(length(y));
  ll = gaussian_prob(y, mu, Sigma, 1);
  %ll = log(gaussian_prob(y, mu, Sigma));
  return
end

[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, dom, sizes);
K = SXY*inv(SYY);
mu2 = muX + K*(y-muY);
Sigma2 = SXX - K*SYX;
ll = gaussian_prob(y, muY, SYY, 1);
%ll = log(gaussian_prob(y, muY, SYY));



