function bnode = maximize_params_gaussian(bnode)
% MAXIMIZE_PARAMS_GAUSSIAN Set the parameters of a Gaussian node to their ML value.
% bnode = maximize_params_gaussian(bnode)
%
% For details, see "Fitting a Conditional Gaussian Distribution", Kevin Murphy, tech. note, 1998.

%fprintf('updating node %d\n', bnode.self);

ss = bnode.self_size;
dpsize = bnode.dpsize;
cpsize = bnode.cpsize;
dsizes = bnode.dpsizes;

check_covars = 1;
check_real = 1;
min_covar = 1e-5;

w = bnode.Wsum(:);
% Set any zeros to one before dividing
% This is valid because w(i)=0 iff WYsum(:,i)=0, etc.
w = w + (w==0);

if check_real
  % Sometimes, with limited data, we can get imaginary components!
  % These are basically harmless, but ugly.
  if ~isreal(bnode.WXYsum) | ~isreal(bnode.WXXsum) | ~isreal(bnode.WYYsum)
    disp(['removing imaginary components for node ' num2str(bnode.self)]);
    bnode.WXYsum = real(bnode.WXYsum);
    bnode.WXXsum = real(bnode.WXXsum);
    bnode.WYYsum = real(bnode.WYYsum);
  end
end

% Let x2 = [x 1]' and B2 = [B mu]
B2 = zeros(ss, cpsize+1, dpsize);
XY = zeros(cpsize+1, ss, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)' 
XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)' 
YY = zeros(ss, ss, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' 
for i=1:dpsize
  XY(:,:,i) = [bnode.WXYsum(:,:,i)
               bnode.WYsum(:,i)'];
  XX(:,:,i) = [bnode.WXXsum(:,:,i) bnode.WXsum(:,i);
               bnode.WXsum(:,i)'   bnode.Wsum(i)];
  YY(:,:,i) = bnode.WYYsum(:,:,i);
end

zero_mean =  ~bnode.adjustable_mean & isequal(bnode.mean, zeros(bnode.self_size, 1, bnode.dpsize));
% Estimating B2 and then setting the last column (the mean) to 0s is *not* equivalent
% to estimating B and then adding 0s to the last column.
if zero_mean
  if bnode.adjustable_weights & cpsize>0 
    B = zeros(ss, cpsize, dpsize);
    for i=1:dpsize
      if det(bnode.WXXsum(:,:,i))==0
	B(:,:,i) = 0;
      else
	%B(:,:,i) = bnode.WXYsum(:,:,i)' * inv(bnode.WXXsum(:,:,i));
	temp = bnode.WXXsum(:,:,i) \ bnode.WXYsum(:,:,i);
	B(:,:,i) = temp';
      end
    end
    bnode.weights = reshape(B, [ss cpsize dsizes]);
  end
else
  % Estimate mean and weights simultaneously
  for i=1:dpsize
    if det(XX(:,:,i))==0  % fix by U. Sondhauss 6/27/99
      B2(:,:,i)=0;          
    else                    
      %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
      temp = XX(:,:,i) \ XY(:,:,i);
      B2(:,:,i) = temp';
    end                   
  end
  for i=1:dpsize
    if bnode.adjustable_mean
      bnode.mean(:,i) = B2(:,cpsize+1,i);
    end
    if bnode.adjustable_weights
      bnode.weights(:,:,i) = B2(:,1:cpsize,i);
    end
  end
end

if cpsize>0
  B2(:,1:cpsize,:) = reshape(bnode.weights, [ss cpsize dpsize]);
end
B2(:,cpsize+1,:) = reshape(bnode.mean, [ss dpsize]);

if bnode.adjustable_cov
  if bnode.tied_cov
    S = zeros(ss, ss);
    for i=1:dpsize
      S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
    end
    S = S / bnode.num_samples;
    if strcmp(bnode.cov_type, 'diag')
      S = diag(diag(S));
    end
    bnode.cov = repmat(S, [1 1 dpsize]);
  else 
    for i=1:dpsize      
      bnode.cov(:,:,i) = (YY(:,:,i) - B2(:,:,i)*XY(:,:,i))/w(i);
    end
    if strcmp(bnode.cov_type, 'diag')
      for i=1:dpsize      
	bnode.cov(:,:,i) = diag(diag(bnode.cov(:,:,i)));
      end
    end
    
  end
end

if check_covars % prevent collapsing to a point
  for i=1:dpsize
    if min(svd(bnode.cov(:,:,i))) < min_covar
      disp(['resetting singular covariance for node ' num2str(bnode.self)]);
      bnode.cov(:,:,i) = bnode.init_cov(:,:,i);
    end
  end
end



