function bnode = mk_gaussian_node(bnet, self, cov_type, tied_cov)
% MK_GAUSSIAN_NODE Make a node whose cond. distrib. is a Gaussian with random params.
% bnode = mk_gaussian_node(bnet, node_num, cov_type, tied_cov)
%
% The type of distribution depends on the kind of parents this node (call it Y) has.
% - no parents: Y ~ N(mu, Sigma)
% - continuous parents X: Y|X=x ~ N(mu + W x, Sigma)
% - discrete parents Q: Y|Q=i ~ N(mu(i), Sigma(i))
% - cts parents X and discrete parents Q: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))
%
% The params are called 'mean(:,i)', 'cov(:,:,i)' and 'weights(:.:.i)' for each value i
% of the discrete parents.
%
% You can control the num. params in Sigma as follows.
% cov_type = 'full' means Sigma(i) is unconstrained 
% cov_type = 'diag' means Sigma(i) is diagonal. Default = 'full'.
% tied_cov = 1 means Sigma(i) = Sigma for all values i of Q. Default = 0.
%
% If bnet.node_sizes(Y)=1, Y will be a scalar, otherwise a vector.
%
% Use randn('state', seed) before calling to control the random param. values.

if nargin < 3, cov_type = 'full'; end
if nargin < 4, tied_cov = 0; end

if ~ismember(self, bnet.cts_nodes)
  error('gaussian nodes must be continuous');
end

bnode.type = 'gaussian';
bnode.self = self;


bnode.parents = parents(bnet.adj_mat, self);
ns = bnet.node_sizes;
ss = ns(self);
bnode.self_size = ss;
bnode.dparents = intersect(bnode.parents, bnet.discrete_nodes);
bnode.cparents = intersect(bnode.parents, bnet.cts_nodes);
bnode.dpsize = prod(ns(bnode.dparents));
bnode.dpsizes = ns(bnode.dparents);
bnode.cpsize = sum(ns(bnode.cparents));
bnode.cfamsize = ns([bnode.cparents bnode.self]);

bnode.adjustable_mean = 1;
bnode.adjustable_cov = 1;
bnode.adjustable_weights = 1;
bnode.adjustable = 1;
bnode.cov_type = cov_type;
bnode.tied_cov = tied_cov;

switch cov_type
  case 'full',
    ncov_params = ss*(ss-1)/2; % since symmetric (and positive definite)
  case 'diag',
    ncov_params = ss;
  otherwise
    error(['unrecognized cov_type ' cov_type]);
end
% params = weights + mean + cov
if bnode.tied_cov
  bnode.nparams = ss*bnode.cpsize*bnode.dpsize + ss*bnode.dpsize + ncov_params;
else
  bnode.nparams = ss*bnode.cpsize*bnode.dpsize + ss*bnode.dpsize + bnode.dpsize*ncov_params;
end

bnode = randomize_params_gaussian(bnode);
bnode = reset_ess_gaussian(bnode);

bnode.adjustable_weights = (bnode.cpsize > 0);
