function bnode = update_ess_gaussian(bnode, marginal)
% UPDATE_ESS_GAUSSIAN Update the Expected Sufficient Statistics for a Gaussian node.
% bnode = update_ess_gaussian(bnode, marginal)

adjustable = bnode.adjustable_mean | bnode.adjustable_cov | bnode.adjustable_weights;
if ~adjustable, return; end

self = bnode.self;

bnode.num_samples = bnode.num_samples + 1;

if isempty(bnode.dparents)
  w = 1;
else
  w = marginal.T(:);
end
bnode.Wsum = bnode.Wsum + w;
% Let X be the cts parent, Y be the cts child (self).
% E[X Y] = Cov[X,Y] + E[X] E[Y]
for i=1:bnode.dpsize
  [muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(marginal.mu(:,i), marginal.Sigma(:,:,i), ...
      bnode.cparents, bnode.self, [bnode.cparents bnode.self], bnode.cfamsize);
  bnode.WYsum(:,i) = bnode.WYsum(:,i) + w(i)*muY;
  bnode.WYYsum(:,:,i) = bnode.WYYsum(:,:,i) + w(i)*(SYY + muY*muY');
  if bnode.cpsize > 0
    bnode.WXsum(:,i) = bnode.WXsum(:,i) + w(i)*muX;
    bnode.WXYsum(:,:,i) = bnode.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
    bnode.WXXsum(:,:,i) = bnode.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');  
  end
end
