BNW using Octave instead of Matlab.

This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning.

I am calling this BNW_1.02. It can be accessed at:
compbio.uthsc.edu/BNW_1.02

36 files changed, 1699 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_lambda_msg.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_lambda_msg.m
new file mode 100644
index 00000000..340ebe5c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_lambda_msg.m
@@ -0,0 +1,59 @@
+function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
+% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)
+% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
+% Pearl p183 eq 4.52
+
+switch msg_type
+ case 'd',
+  error('gaussian_CPD can''t create discrete msgs')
+ case 'g',
+  cps = ps(CPD.cps);
+  cpsizes = CPD.sizes(CPD.cps);
+  self_size = CPD.sizes(end);
+  i = find_equiv_posns(p, cps); % p is n's i'th cts parent
+  psz = cpsizes(i);
+  if all(msg{n}.lambda.precision == 0) % no info to send on
+    lam_msg.precision = zeros(psz, psz);
+    lam_msg.info_state = zeros(psz, 1);
+    return;
+  end
+  [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence);
+  Bmu = m;
+  BSigma = Q;
+  for k=1:length(cps) % only get pi msgs from cts parents
+    pk = cps(k);
+    if pk ~= p
+      %bk = block(k, cpsizes);
+      bk = CPD.cps_block_ndx{k};
+      Bk = W(:, bk);
+      m = msg{n}.pi_from_parent{k}; 
+      BSigma = BSigma + Bk * m.Sigma * Bk';
+      Bmu = Bmu + Bk * m.mu;
+    end
+  end
+  % BSigma = Q + sum_{k \neq i} B_k Sigma_k B_k'
+  %bi = block(i, cpsizes);
+  bi = CPD.cps_block_ndx{i};
+  Bi = W(:,bi);
+  P = msg{n}.lambda.precision;
+  if (rcond(P) > 1e-3) || isinf(P)
+    if isinf(P) % Y is observed
+      Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance
+      mu_lambda = msg{n}.lambda.mu; % observed_value;
+    else
+      Sigma_lambda = inv(P);
+      mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;
+    end
+    C = inv(Sigma_lambda + BSigma);
+    lam_msg.precision = Bi' * C * Bi;
+    lam_msg.info_state = Bi' * C * (mu_lambda - Bmu);
+  else
+    % method that uses matrix inversion lemma to avoid inverting P
+    A = inv(P + inv(BSigma));
+    C = P - P*A*P;
+    lam_msg.precision = Bi' * C * Bi;
+    D = eye(self_size) - P*A;
+    z = msg{n}.lambda.info_state;
+    lam_msg.info_state = Bi' * (D*z - D*P*Bmu);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_pi.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_pi.m
new file mode 100644
index 00000000..910973e7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_pi.m
@@ -0,0 +1,22 @@
+function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)
+% CPD_TO_PI Compute the pi vector (gaussian)
+% function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)
+
+switch msg_type
+ case 'd',
+  error('gaussian_CPD can''t create discrete msgs')
+ case 'g',
+  [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence);
+  cps = ps(CPD.cps);
+  cpsizes = CPD.sizes(CPD.cps);
+  pi.mu = m;
+  pi.Sigma = Q;
+  for k=1:length(cps) % only get pi msgs from cts parents
+    %bk = block(k, cpsizes);
+    bk = CPD.cps_block_ndx{k};
+    Bk = W(:, bk);
+    m = msg{n}.pi_from_parent{k}; 
+    pi.Sigma = pi.Sigma + Bk * m.Sigma * Bk';
+    pi.mu = pi.mu + Bk * m.mu; % m.mu = u(k)
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_scgpot.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_scgpot.m
new file mode 100644
index 00000000..90e7cc80
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CPD_to_scgpot.m
@@ -0,0 +1,58 @@
+function pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence)
+% CPD_TO_CGPOT Convert a Gaussian CPD to a CG potential, incorporating any evidence   
+% pot = CPD_to_cgpot(CPD, domain, ns, cnodes, evidence)
+
+self = CPD.self;
+dnodes = mysetdiff(1:length(ns), cnodes);
+odom = domain(~isemptycell(evidence(domain)));
+cdom = myintersect(cnodes, domain);
+cheaddom = myintersect(self, domain);
+ctaildom = mysetdiff(cdom,cheaddom);
+ddom = myintersect(dnodes, domain);
+cobs = myintersect(cdom, odom);
+dobs = myintersect(ddom, odom);
+ens = ns; % effective node size
+ens(cobs) = 0;
+ens(dobs) = 1;
+
+% Extract the params compatible with the observations (if any) on the discrete parents (if any)
+% parents are all but the last domain element
+ps = domain(1:end-1);
+dps = myintersect(ps, ddom);
+dops = myintersect(dps, odom);
+
+map = find_equiv_posns(dops, dps);
+dpvals = cat(1, evidence{dops});
+index = mk_multi_index(length(dps), map, dpvals);
+
+dpsize = prod(ens(dps));
+cpsize = size(CPD.weights(:,:,1), 2); % cts parents size
+ss = size(CPD.mean, 1); % self size
+% the reshape acts like a squeeze
+m = reshape(CPD.mean(:, index{:}), [ss dpsize]);
+C = reshape(CPD.cov(:, :, index{:}), [ss ss dpsize]);
+W = reshape(CPD.weights(:, :, index{:}), [ss cpsize dpsize]);
+
+
+% Convert each conditional Gaussian to a canonical potential
+pot = cell(1, dpsize);
+for i=1:dpsize
+  %pot{i} = linear_gaussian_to_scgcpot(m(:,i), C(:,:,i), W(:,:,i), cdom, ns, cnodes, evidence);
+  pot{i} = scgcpot(ss, cpsize, 1, m(:,i), W(:,:,i), C(:,:,i));
+end
+
+pot = scgpot(ddom, cheaddom, ctaildom, ens, pot);
+
+
+function pot = linear_gaussian_to_scgcpot(mu, Sigma, W, domain, ns, cnodes, evidence)
+% LINEAR_GAUSSIAN_TO_CPOT Convert a linear Gaussian CPD  to a stable conditional potential element.
+% pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence)
+
+p = 1;
+A = mu;
+B = W;
+C = Sigma;
+ns(odom) = 0;
+%pot = scgcpot(, ns(domain), p, A, B, C);
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Entries
new file mode 100644
index 00000000..a6bd3e14
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Entries
@@ -0,0 +1,20 @@
+/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/CPD_to_scgpot.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/adjustable_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_CPD_to_table_hidden_ps.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/convert_to_pot.m/1.1.1.1/Sun Mar  9 23:03:16 2003//
+/convert_to_table.m/1.1.1.1/Sun May 11 23:31:54 2003//
+/display.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/gaussian_CPD.m/1.1.1.1/Wed Jun 15 21:13:06 2005//
+/gaussian_CPD_params_given_dps.m/1.1.1.1/Sun May 11 23:13:40 2003//
+/get_field.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/learn_params.m/1.1.1.1/Thu Jun 10 01:28:10 2004//
+/log_prob_node.m/1.1.1.1/Tue Sep 10 17:44:00 2002//
+/maximize_params.m/1.1.1.1/Tue May 20 14:10:06 2003//
+/maximize_params_debug.m/1.1.1.1/Fri Jan 31 00:13:10 2003//
+/reset_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/sample_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/set_fields.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/update_ess.m/1.1.1.1/Tue Jul 22 22:55:46 2003//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Entries.Log
new file mode 100644
index 00000000..9c6f22e4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Entries.Log
@@ -0,0 +1,2 @@
+A D/Old////
+A D/private////
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Repository
new file mode 100644
index 00000000..98ebf3cb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@gaussian_CPD
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m
new file mode 100644
index 00000000..5a6d398a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m
@@ -0,0 +1,64 @@
+function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p)
+% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)
+% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p)
+% Pearl p183 eq 4.52
+
+switch msg_type
+ case 'd',
+  error('gaussian_CPD can''t create discrete msgs')
+ case 'g',
+  self_size = CPD.sizes(end);
+  if all(msg{n}.lambda.precision == 0) % no info to send on
+    lam_msg.precision = zeros(self_size);
+    lam_msg.info_state = zeros(self_size, 1);
+    return;
+  end
+  cpsizes = CPD.sizes(CPD.cps);
+  dpval = 1;
+  Q = CPD.cov(:,:,dpval);
+  Sigmai = Q;
+  wmu = zeros(self_size, 1);
+  for k=1:length(ps)
+    pk = ps(k);
+    if pk ~= p
+      bk = block(k, cpsizes);
+      Bk = CPD.weights(:, bk, dpval);
+      m = msg{n}.pi_from_parent{k};
+      Sigmai = Sigmai + Bk * m.Sigma * Bk';
+      wmu = wmu + Bk * m.mu; % m.mu = u(k)
+    end
+  end
+  % Sigmai = Q + sum_{k \neq i} B_k Sigma_k B_k'
+  i = find_equiv_posns(p, ps);
+  bi = block(i, cpsizes);
+  Bi = CPD.weights(:,bi, dpval);
+  
+  if 0
+  P = msg{n}.lambda.precision;
+  if isinf(P) % inv(P)=Sigma_lambda=0
+    precision_temp = inv(Sigmai);
+    lam_msg.precision = Bi' * precision_temp * Bi;
+    lam_msg.info_state = precision_temp * (msg{n}.lambda.mu - wmu);
+  else
+    A = inv(P + inv(Sigmai));
+    precision_temp = P + P*A*P;
+    lam_msg.precision = Bi' * precision_temp * Bi;
+    self_size = length(P);
+    C = eye(self_size) + P*A;
+    z = msg{n}.lambda.info_state;
+    lam_msg.info_state = C*z - C*P*wmu;
+  end
+  end
+  
+  if isinf(msg{n}.lambda.precision)
+    Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance
+    mu_lambda = msg{n}.lambda.mu; % observed_value;
+  else
+    Sigma_lambda = inv(msg{n}.lambda.precision);
+    mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;
+  end
+  precision_temp = inv(Sigma_lambda + Sigmai);
+  lam_msg.precision = Bi' * precision_temp * Bi;
+  lam_msg.info_state = Bi' * precision_temp * (mu_lambda - wmu);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries
new file mode 100644
index 00000000..ea2f5a4c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries
@@ -0,0 +1,7 @@
+/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/gaussian_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/maximize_params.m/1.1.1.1/Thu Jan 30 22:38:16 2003//
+/update_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/update_tied_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository
new file mode 100644
index 00000000..c89b5b86
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@gaussian_CPD/Old
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m
new file mode 100644
index 00000000..6f7138fc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m
@@ -0,0 +1,184 @@
+function CPD = gaussian_CPD(varargin)
+% GAUSSIAN_CPD Make a conditional linear Gaussian distrib.
+%
+% To define this CPD precisely, call the continuous (cts) parents (if any) X,
+% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is:
+% - no parents: Y ~ N(mu, Sigma)
+% - cts parents : Y|X=x ~ N(mu + W x, Sigma)
+% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i))
+% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))
+%
+% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters,
+% where node is the number of a node in this equivalence class.
+%
+% The list below gives optional arguments [default value in brackets].
+% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).)
+%
+% mean       - mu(:,i) is the mean given Q=i [ randn(Y,Q) ]
+% cov        - Sigma(:,:,i) is the covariance given Q=i [ repmat(eye(Y,Y), [1 1 Q]) ]
+% weights    - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ]
+% cov_type   - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ]
+% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i [0]
+% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0]
+% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning [0]
+% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0]
+% cov_prior_weight - weight given to I prior for estimating Sigma [0.01]
+%
+% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 'yes')
+%
+% For backwards compatibility with BNT2, you can also specify the parameters in the following order
+%   CPD = gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weight)
+%
+% Sometimes it is useful to create an "isolated" CPD, without needing to pass in a bnet.
+% In this case, you must specify the discrete and cts parents (dps, cps) and the family sizes, followed
+% by the optional arguments above:
+%   CPD = gaussian_CPD('self', i, 'dps', dps, 'cps', cps, 'sz', fam_size, ...)
+
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  clamp = 0;
+  CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp));
+  return;
+elseif isa(varargin{1}, 'gaussian_CPD')
+  % This might occur if we are copying an object.
+  CPD = varargin{1};
+  return;
+end
+CPD = init_fields;
+ 
+CPD = class(CPD, 'gaussian_CPD', generic_CPD(0));
+
+
+% parse mandatory arguments
+if ~isstr(varargin{1}) % pass in bnet
+  bnet = varargin{1};
+  self = varargin{2};
+  args = varargin(3:end);
+  ns = bnet.node_sizes;
+  ps = parents(bnet.dag, self);
+  dps = myintersect(ps, bnet.dnodes);
+  cps = myintersect(ps, bnet.cnodes);
+  fam_sz = ns([ps self]);
+else
+  disp('parsing new style')
+  for i=1:2:length(varargin)
+    switch varargin{i},
+     case 'self', self = varargin{i+1}; 
+     case 'dps',  dps = varargin{i+1};
+     case 'cps',  cps = varargin{i+1};
+     case 'sz',   fam_sz = varargin{i+1};
+    end
+  end
+  ps = myunion(dps, cps);
+  args = varargin;
+end
+
+CPD.self = self;
+CPD.sizes = fam_sz;
+
+% Figure out which (if any) of the parents are discrete, and which cts, and how big they are
+% dps = discrete parents, cps = cts parents
+CPD.cps = find_equiv_posns(cps, ps); % cts parent index
+CPD.dps = find_equiv_posns(dps, ps);
+ss = fam_sz(end);
+psz = fam_sz(1:end-1);
+dpsz = prod(psz(CPD.dps));
+cpsz = sum(psz(CPD.cps));
+
+% set default params
+CPD.mean = randn(ss, dpsz);
+CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]);    
+CPD.weights = randn(ss, cpsz, dpsz);
+CPD.cov_type = 'full';
+CPD.tied_cov = 0;
+CPD.clamped_mean = 0;
+CPD.clamped_cov = 0;
+CPD.clamped_weights = 0;
+CPD.cov_prior_weight = 0.01;
+
+nargs = length(args);
+if nargs > 0
+  if ~isstr(args{1})
+    % gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weights)
+    if nargs >= 1 & ~isempty(args{1}), CPD.mean = args{1}; end
+    if nargs >= 2 & ~isempty(args{2}), CPD.cov = args{2}; end
+    if nargs >= 3 & ~isempty(args{3}), CPD.weights = args{3}; end
+    if nargs >= 4 & ~isempty(args{4}), CPD.cov_type = args{4}; end
+    if nargs >= 5 & ~isempty(args{5}) & strcmp(args{5}, 'tied'), CPD.tied_cov = 1; end
+    if nargs >= 6 & ~isempty(args{6}), CPD.clamped_mean = 1; end
+    if nargs >= 7 & ~isempty(args{7}), CPD.clamped_cov = 1; end
+    if nargs >= 8 & ~isempty(args{8}), CPD.clamped_weights = 1; end
+  else
+    CPD = set_fields(CPD, args{:});
+  end
+end
+
+% Make sure the matrices have 1 dimension per discrete parent.
+% Bug fix due to Xuejing Sun 3/6/01
+CPD.mean = myreshape(CPD.mean, [ss ns(dps)]);
+CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]);
+CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]);
+  
+CPD.init_cov = CPD.cov;  % we reset to this if things go wrong during learning
+
+% expected sufficient statistics 
+CPD.Wsum = zeros(dpsz,1);
+CPD.WYsum = zeros(ss, dpsz);
+CPD.WXsum = zeros(cpsz, dpsz);
+CPD.WYYsum = zeros(ss, ss, dpsz);
+CPD.WXXsum = zeros(cpsz, cpsz, dpsz);
+CPD.WXYsum = zeros(cpsz, ss, dpsz);
+
+% For BIC
+CPD.nsamples = 0;
+switch CPD.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 CPD.tied_cov
+  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params;
+else
+  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params;
+end
+
+
+
+clamped = CPD.clamped_mean & CPD.clamped_cov & CPD.clamped_weights;
+CPD = set_clamped(CPD, clamped);
+
+%%%%%%%%%%%
+
+function CPD = init_fields()
+% This ensures we define the fields in the same order 
+% no matter whether we load an object from a file,
+% or create it from scratch. (Matlab requires this.)
+
+CPD.self = [];
+CPD.sizes = [];
+CPD.cps = [];
+CPD.dps = [];
+CPD.mean = [];
+CPD.cov = [];
+CPD.weights = [];
+CPD.clamped_mean = [];
+CPD.clamped_cov = [];
+CPD.clamped_weights = [];
+CPD.init_cov = [];
+CPD.cov_type = [];
+CPD.tied_cov = [];
+CPD.Wsum = [];
+CPD.WYsum = [];
+CPD.WXsum = [];
+CPD.WYYsum = [];
+CPD.WXXsum = [];
+CPD.WXYsum = [];
+CPD.nsamples = [];
+CPD.nparams = [];            
+CPD.cov_prior_weight = [];
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m
new file mode 100644
index 00000000..3fa398c8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m
@@ -0,0 +1,59 @@
+function L = log_prob_node(CPD, self_ev, pev)
+% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (gaussian)
+% L = log_prob_node(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+
+if iscell(self_ev), usecell = 1; else usecell = 0; end
+
+use_log = 1;
+ncases = length(self_ev);
+nparents = length(CPD.sizes)-1;
+assert(ncases == size(pev, 2));
+
+if ncases == 0
+  L = 0;
+  return;
+end
+
+if length(CPD.dps)==0 % no discrete parents, so we can vectorize
+  i = 1;
+  if usecell
+    Y = cell2num(self_ev);
+  else
+    Y = self_ev;
+  end
+  if length(CPD.cps) == 0 
+    L = gaussian_prob(Y, CPD.mean(:,i), CPD.cov(:,:,i), use_log);
+  else
+    if usecell
+      X = cell2num(pev);
+    else
+      X = pev;
+    end
+    L = gaussian_prob(Y, CPD.mean(:,i) + CPD.weights(:,:,i)*X, CPD.cov(:,:,i), use_log);
+  end
+else % each case uses a (potentially) different set of parameters
+  L = 0;
+  for m=1:ncases
+    if usecell
+      dpvals = cat(1, pev{CPD.dps, m});
+    else
+      dpvals = pev(CPD.dps, m);
+    end
+    i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)');
+    y = self_ev{m};
+    if length(CPD.cps) == 0 
+      L = L + gaussian_prob(y, CPD.mean(:,i), CPD.cov(:,:,i), use_log);
+    else
+      if usecell
+	x = cat(1, pev{CPD.cps, m});
+      else
+	x = pev(CPD.cps, m);
+      end
+      L = L + gaussian_prob(y, CPD.mean(:,i) + CPD.weights(:,:,i)*x, CPD.cov(:,:,i), use_log);
+    end
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m
new file mode 100644
index 00000000..48447358
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m
@@ -0,0 +1,147 @@
+function CPD = maximize_params(CPD, temp)
+% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
+% CPD = maximize_params(CPD, temperature)
+%
+% Temperature is currently only used for entropic prior on Sigma
+
+% For details, see "Fitting a Conditional Gaussian Distribution", Kevin Murphy, tech. report,
+% 1998, available at www.cs.berkeley.edu/~murphyk/papers.html
+% Refering to table 2, we use equations 1/2 to estimate the covariance matrix in the untied/tied case,
+% and equation 9 to estimate the weight matrix and mean.
+% We do not implement spherical Gaussians - the code is already pretty complicated!
+
+if ~adjustable_CPD(CPD), return; end
+
+%assert(approxeq(CPD.nsamples, sum(CPD.Wsum)));
+assert(~any(isnan(CPD.WXXsum)))
+assert(~any(isnan(CPD.WXYsum)))
+assert(~any(isnan(CPD.WYYsum)))
+
+[self_size cpsize dpsize] = size(CPD.weights);
+
+% Append 1s to the parents, and derive the corresponding cross products.
+% This is used when estimate the means and weights simultaneosuly,
+% and when estimatting Sigma.
+% Let x2 = [x 1]'
+XY = zeros(cpsize+1, self_size, 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(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' 
+for i=1:dpsize
+  XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y
+	       CPD.WYsum(:,i)']; % 1*Y
+  % [x  * [x' 1]  = [xx' x
+  %  1]              x'  1]
+  XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i);
+	       CPD.WXsum(:,i)'   CPD.Wsum(i)];
+  YY(:,:,i) = CPD.WYYsum(:,:,i);
+end
+
+w = CPD.Wsum(:);
+% Set any zeros to one before dividing
+% This is valid because w(i)=0 => WYsum(:,i)=0, etc
+w = w + (w==0);
+
+if CPD.clamped_mean
+  % Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent
+  % to estimating B and then adding the clamped_mean to the last column.
+  if ~CPD.clamped_weights
+    B = zeros(self_size, cpsize, dpsize);
+    for i=1:dpsize
+      if det(CPD.WXXsum(:,:,i))==0
+	B(:,:,i) = 0;
+      else
+	% Eqn 9 in table 2 of TR
+	%B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i));
+	B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))';
+      end
+    end
+    %CPD.weights = reshape(B, [self_size cpsize dpsize]);
+    CPD.weights = B;
+  end
+elseif CPD.clamped_weights % KPM 1/25/02
+  if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals
+    for i=1:dpsize
+      CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i);
+    end
+  end
+else % nothing is clamped, so estimate mean and weights simultaneously
+  B2 = zeros(self_size, cpsize+1, dpsize);
+  for i=1:dpsize
+    if det(XX(:,:,i))==0  % fix by U. Sondhauss 6/27/99
+      B2(:,:,i)=0;          
+    else                    
+      % Eqn 9 in table 2 of TR
+      %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
+      B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))';
+    end                   
+    CPD.mean(:,i) = B2(:,cpsize+1,i);
+    CPD.weights(:,:,i) = B2(:,1:cpsize,i);
+  end
+end
+
+% Let B2 = [W mu]
+if cpsize>0
+  B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]);
+end
+B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]);
+
+% To avoid singular covariance matrices,
+% we use the regularization method suggested in "A Quasi-Bayesian approach to estimating
+% parameters for mixtures of normal distributions", Hamilton 91.
+% If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma),
+% where gamma >=0 is a smoothing parameter (equivalent sample size of I prior)
+
+gamma = CPD.cov_prior_weight;
+
+if ~CPD.clamped_cov
+  if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99
+    Z = 1-temp;
+    % When temp > 1, Z is negative, so we are dividing by a smaller
+    % number, ie. increasing the variance.
+  else
+    Z = 0;
+  end
+  if CPD.tied_cov
+    S = zeros(self_size, self_size);
+    % Eqn 2 from table 2 in TR
+    for i=1:dpsize
+      S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
+    end
+    %denom = max(1, CPD.nsamples + gamma + Z);
+    denom = CPD.nsamples + gamma + Z;
+    S = (S + gamma*eye(self_size)) / denom;
+    if strcmp(CPD.cov_type, 'diag')
+      S = diag(diag(S));
+    end
+    CPD.cov = repmat(S, [1 1 dpsize]);
+  else 
+    for i=1:dpsize      
+      % Eqn 1 from table 2 in TR
+      S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i);
+      %denom = max(1, w(i) + gamma + Z); % gives wrong answers on mhmm1
+      denom = w(i) + gamma + Z;
+      S = (S + gamma*eye(self_size)) / denom;
+      CPD.cov(:,:,i) = S;
+    end
+    if strcmp(CPD.cov_type, 'diag')
+      for i=1:dpsize      
+	CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i)));
+      end
+    end
+  end
+end
+
+
+check_covars = 0;
+min_covar = 1e-5;
+if check_covars % prevent collapsing to a point
+  for i=1:dpsize
+    if min(svd(CPD.cov(:,:,i))) < min_covar
+      disp(['resetting singular covariance for node ' num2str(CPD.self)]);
+      CPD.cov(:,:,i) = CPD.init_cov(:,:,i);
+    end
+  end
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_ess.m
new file mode 100644
index 00000000..988012e2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_ess.m
@@ -0,0 +1,85 @@
+function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+% UPDATE_ESS Update the Expected Sufficient Statistics of a Gaussian node
+% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+
+%if nargin < 6
+%  hidden_bitv = zeros(1, max(fmarginal.domain));
+%  hidden_bitv(find(isempty(evidence)))=1;
+%end
+
+dom = fmarginal.domain;
+self = dom(end);
+ps = dom(1:end-1);
+hidden_self = hidden_bitv(self);
+cps = myintersect(ps, cnodes);
+dps = mysetdiff(ps, cps);
+hidden_cps = all(hidden_bitv(cps));
+hidden_dps = all(hidden_bitv(dps));
+
+CPD.nsamples = CPD.nsamples + 1;            
+[ss cpsz dpsz] = size(CPD.weights); % ss = self size
+
+% Let X be the cts parent (if any), Y be the cts child (self).
+
+if ~hidden_self & (isempty(cps) | ~hidden_cps) & hidden_dps % all cts nodes are observed, all discrete nodes are hidden
+  % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0
+  % Since discrete parents are hidden, we do not need to add evidence to w.
+  w = fmarginal.T(:);
+  CPD.Wsum = CPD.Wsum + w;
+  y = evidence{self};
+  Cyy = y*y';
+  if ~CPD.useC
+     W = repmat(w(:)',ss,1); % W(y,i) = w(i)
+     W2 = repmat(reshape(W, [ss 1 dpsz]), [1 ss 1]); % W2(x,y,i) = w(i)
+     CPD.WYsum = CPD.WYsum +  W .* repmat(y(:), 1, dpsz);
+     CPD.WYYsum = CPD.WYYsum + W2  .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]);
+  else
+     W = w(:)';
+     W2 = reshape(W, [1 1 dpsz]);
+     CPD.WYsum = CPD.WYsum +  rep_mult(W, y(:), size(CPD.WYsum)); 
+     CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum));
+  end
+  if cpsz > 0 % X exists
+    x = cat(1, evidence{cps}); x = x(:);
+    Cxx = x*x';
+    Cxy = x*y';
+    if ~CPD.useC
+       CPD.WXsum = CPD.WXsum + W .* repmat(x(:), 1, dpsz);
+       CPD.WXXsum = CPD.WXXsum + W2 .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]);
+       CPD.WXYsum = CPD.WXYsum + W2 .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]);
+    else
+       CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum));
+       CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum));
+       CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum));
+    end
+  end
+  return;
+end
+
+% general (non-vectorized) case
+fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow!
+
+if dpsz == 1 % no discrete parents
+  w = 1;
+else
+  w = fullm.T(:);
+end
+
+CPD.Wsum = CPD.Wsum + w;
+xi = 1:cpsz;
+yi = (cpsz+1):(cpsz+ss);
+for i=1:dpsz
+  muY = fullm.mu(yi, i);
+  SYY = fullm.Sigma(yi, yi, i);
+  CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
+  CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
+  if cpsz > 0
+    muX = fullm.mu(xi, i);
+    SXX = fullm.Sigma(xi, xi, i);
+    SXY = fullm.Sigma(xi, yi, i);
+    CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
+    CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
+    CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
+  end
+end                
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m
new file mode 100644
index 00000000..798c795c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m
@@ -0,0 +1,118 @@
+function CPD = update_tied_ess(CPD, domain, engine, evidence, ns, cnodes)
+
+if ~adjustable_CPD(CPD), return; end
+nCPDs = size(domain, 2);
+fmarginal = cell(1, nCPDs);
+for l=1:nCPDs
+  fmarginal{l} = marginal_family(engine, nodes(l));
+end
+
+[ss cpsz dpsz] = size(CPD.weights);
+if const_evidence_pattern(engine)
+  dom = domain(:,1);
+  dnodes = mysetdiff(1:length(ns), cnodes);
+  ddom = myintersect(dom, dnodes);
+  cdom = myintersect(dom, cnodes);
+  odom = dom(~isemptycell(evidence(dom)));
+  hdom = dom(isemptycell(evidence(dom)));
+  % If all hidden nodes are discrete and all cts nodes are observed 
+  % (e.g., HMM with Gaussian output)
+  % we can add the observed evidence in parallel
+  if mysubset(ddom, hdom) & mysubset(cdom, odom)
+    [mu, Sigma, T] = add_cts_ev_to_marginals(fmarginal, evidence, ns, cnodes);
+  else
+    mu = zeros(ss, dpsz, nCPDs);
+    Sigma = zeros(ss, ss, dpsz, nCPDs);
+    T = zeros(dpsz, nCPDs);
+    for l=1:nCPDs
+      [mu(:,:,l), Sigma(:,:,:,l), T(:,l)] = add_ev_to_marginals(fmarginal{l}, evidence, ns, cnodes);
+    end
+  end
+end
+CPD.nsamples = CPD.nsamples + nCPDs;            
+
+
+if dpsz == 1 % no discrete parents
+  w = 1;
+else
+  w = fullm.T(:);
+end
+CPD.Wsum = CPD.Wsum + w;
+% Let X be the cts parent (if any), Y be the cts child (self).
+xi = 1:cpsz;
+yi = (cpsz+1):(cpsz+ss);
+for i=1:dpsz
+  muY = fullm.mu(yi, i);
+  SYY = fullm.Sigma(yi, yi, i);
+  CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
+  CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
+  if cpsz > 0
+    muX = fullm.mu(xi, i);
+    SXX = fullm.Sigma(xi, xi, i);
+    SXY = fullm.Sigma(xi, yi, i);
+    CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
+    CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
+    CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
+  end
+end                
+
+
+%%%%%%%%%%%%%
+
+function fullm = add_evidence_to_marginal(fmarginal, evidence, ns, cnodes)
+
+
+dom = fmarginal.domain;
+
+% Find out which values of the discrete parents (if any) are compatible with 
+% the discrete evidence (if any).
+dnodes = mysetdiff(1:length(ns), cnodes);
+ddom = myintersect(dom, dnodes);
+cdom = myintersect(dom, cnodes);
+odom = dom(~isemptycell(evidence(dom)));
+hdom = dom(isemptycell(evidence(dom)));
+
+dobs = myintersect(ddom, odom);
+dvals = cat(1, evidence{dobs});
+ens = ns; % effective node sizes
+ens(dobs) = 1;
+S = prod(ens(ddom));
+subs = ind2subv(ens(ddom), 1:S);
+mask = find_equiv_posns(dobs, ddom);
+subs(mask) = dvals;
+supportedQs = subv2ind(ns(ddom), subs);
+
+if isempty(ddom)
+  Qarity = 1;
+else
+  Qarity = prod(ns(ddom));
+end
+fullm.T = zeros(Qarity, 1);
+fullm.T(supportedQs) = fmarginal.T(:);
+
+% Now put the hidden cts parts into their right blocks,
+% leaving the observed cts parts as 0.
+cobs = myintersect(cdom, odom);
+chid = myintersect(cdom, hdom);
+cvals = cat(1, evidence{cobs});
+n = sum(ns(cdom));
+fullm.mu = zeros(n,Qarity);
+fullm.Sigma = zeros(n,n,Qarity);
+
+if ~isempty(chid)
+  chid_blocks = block(find_equiv_posns(chid, cdom), ns(cdom));
+end
+if ~isempty(cobs)
+  cobs_blocks = block(find_equiv_posns(cobs, cdom), ns(cdom));
+end
+
+for i=1:length(supportedQs)
+  Q = supportedQs(i);
+  if ~isempty(chid)
+    fullm.mu(chid_blocks, Q) = fmarginal.mu(:, i);
+    fullm.Sigma(chid_blocks, chid_blocks, Q) = fmarginal.Sigma(:,:,i);
+  end
+  if ~isempty(cobs)
+    fullm.mu(cobs_blocks, Q) = cvals(:);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/adjustable_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/adjustable_CPD.m
new file mode 100644
index 00000000..ea5190c3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/adjustable_CPD.m
@@ -0,0 +1,5 @@
+function p = adjustable_CPD(CPD)
+% ADJUSTABLE_CPD Does this CPD have any adjustable params? (gaussian)
+% p = adjustable_CPD(CPD)
+
+p = ~CPD.clamped_mean || ~CPD.clamped_cov || ~CPD.clamped_weights;
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_CPD_to_table_hidden_ps.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_CPD_to_table_hidden_ps.m
new file mode 100644
index 00000000..acb2c7d2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_CPD_to_table_hidden_ps.m
@@ -0,0 +1,20 @@
+function T = convert_CPD_to_table_hidden_ps(CPD, self_val)
+% CONVERT_CPD_TO_TABLE_HIDDEN_PS Convert a Gaussian CPD to a table
+% function T = convert_CPD_to_table_hidden_ps(CPD, self_val)
+%
+% self_val must be a non-empty vector.
+% All the parents are hidden.
+%
+% This is used by misc/convert_dbn_CPDs_to_tables
+
+m = CPD.mean;
+C = CPD.cov;
+W = CPD.weights;
+
+[ssz dpsize] = size(m);
+
+T = zeros(dpsize, 1);
+for i=1:dpsize
+  T(i) = gaussian_prob(self_val, m(:,i), C(:,:,i));
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_to_pot.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_to_pot.m
new file mode 100644
index 00000000..6afe8d1d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_to_pot.m
@@ -0,0 +1,71 @@
+function pot = convert_to_pot(CPD, pot_type, domain, evidence)
+% CONVERT_TO_POT Convert a Gaussian CPD to one or more potentials
+% pot = convert_to_pot(CPD, pot_type, domain, evidence)
+
+sz = CPD.sizes;
+ns = zeros(1, max(domain));
+ns(domain) = sz;
+
+odom = domain(~isemptycell(evidence(domain)));
+ps = domain(1:end-1);
+cps = ps(CPD.cps);
+dps = ps(CPD.dps);
+self = domain(end);
+cdom = [cps(:)' self];
+ddom = dps;
+cnodes = cdom;
+  
+switch pot_type
+ case 'u',
+  error('gaussian utility potentials not yet supported');
+ 
+ case 'd',
+  T = convert_to_table(CPD, domain, evidence);
+  ns(odom) = 1;
+  pot = dpot(domain, ns(domain), T);          
+
+ case {'c','g'},
+  [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence);
+  pot = linear_gaussian_to_cpot(m, C, W, domain, ns, cnodes, evidence);
+
+ case 'cg',
+  [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence);
+  % Convert each conditional Gaussian to a canonical potential
+  cobs = myintersect(cdom, odom);
+  dobs = myintersect(ddom, odom);
+  ens = ns; % effective node size
+  ens(cobs) = 0;
+  ens(dobs) = 1;
+  dpsize = prod(ens(dps));
+  can = cell(1, dpsize);
+  for i=1:dpsize
+    if isempty(W)
+      can{i} = linear_gaussian_to_cpot(m(:,i), C(:,:,i), [], cdom, ns, cnodes, evidence);
+    else
+      can{i} = linear_gaussian_to_cpot(m(:,i), C(:,:,i), W(:,:,i), cdom, ns, cnodes, evidence);
+    end
+  end
+  pot = cgpot(ddom, cdom, ens, can);
+
+ case 'scg',
+  [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence);
+  cobs = myintersect(cdom, odom);
+  dobs = myintersect(ddom, odom);
+  ens = ns; % effective node size
+  ens(cobs) = 0;
+  ens(dobs) = 1;
+  dpsize = prod(ens(dps));
+  cpsize = size(W, 2); % cts parents size
+  ss = size(m, 1); % self size
+  cheaddom = self;
+  ctaildom = cps(:)';
+  pot_array = cell(1, dpsize);
+  for i=1:dpsize
+    pot_array{i} = scgcpot(ss, cpsize, 1, m(:,i), W(:,:,i), C(:,:,i));
+  end
+  pot = scgpot(ddom, cheaddom, ctaildom, ens, pot_array);
+
+ otherwise,
+  error(['unrecognized pot_type' pot_type])
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_to_table.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_to_table.m
new file mode 100644
index 00000000..4a8d5904
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/convert_to_table.m
@@ -0,0 +1,38 @@
+function T = convert_to_table(CPD, domain, evidence)
+% CONVERT_TO_TABLE Convert a Gaussian CPD to a table
+% T = convert_to_table(CPD, domain, evidence)
+
+
+sz = CPD.sizes;
+ns = zeros(1, max(domain));
+ns(domain) = sz;
+
+odom = domain(~isemptycell(evidence(domain)));
+ps = domain(1:end-1);
+cps = ps(CPD.cps);
+dps = ps(CPD.dps);
+self = domain(end);
+cdom = [cps(:)' self];
+ddom = dps;
+cnodes = cdom;
+
+[m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence);
+
+
+ns(odom) = 1;
+dpsize = prod(ns(dps));
+self = domain(end);
+assert(myismember(self, odom));
+self_val = evidence{self};
+T = zeros(dpsize, 1);
+if length(cps) > 0 
+  assert(~any(isemptycell(evidence(cps))));
+  cps_vals = cat(1, evidence{cps});
+  for i=1:dpsize
+    T(i) = gaussian_prob(self_val, m(:,i) + W(:,:,i)*cps_vals, C(:,:,i));
+  end
+else
+  for i=1:dpsize
+    T(i) = gaussian_prob(self_val, m(:,i), C(:,:,i));
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/display.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/display.m
new file mode 100644
index 00000000..a3d73c83
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/display.m
@@ -0,0 +1,4 @@
+function display(CPD)
+
+disp('gaussian_CPD object');
+disp(struct(CPD)); 
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/gaussian_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/gaussian_CPD.m
new file mode 100644
index 00000000..de519218
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/gaussian_CPD.m
@@ -0,0 +1,161 @@
+function CPD = gaussian_CPD(bnet, self, varargin)
+% GAUSSIAN_CPD Make a conditional linear Gaussian distrib.
+%
+% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters,
+% where node is the number of a node in this equivalence class.
+
+% To define this CPD precisely, call the continuous (cts) parents (if any) X,
+% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is:
+% - no parents: Y ~ N(mu, Sigma)
+% - cts parents : Y|X=x ~ N(mu + W x, Sigma)
+% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i))
+% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))
+%
+% The list below gives optional arguments [default value in brackets].
+% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).)
+% Parameters will be reshaped to the right size if necessary.
+%
+% mean       - mu(:,i) is the mean given Q=i [ randn(Y,Q) ]
+% cov        - Sigma(:,:,i) is the covariance given Q=i [ repmat(100*eye(Y,Y), [1 1 Q]) ]
+% weights    - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ]
+% cov_type   - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ]
+% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i [0]
+% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0]
+% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning [0]
+% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0]
+% cov_prior_weight - weight given to I prior for estimating Sigma [0.01]
+% cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0]
+%
+% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 1)
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  clamp = 0;
+  CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp));
+  return;
+elseif isa(bnet, 'gaussian_CPD')
+  % This might occur if we are copying an object.
+  CPD = bnet;
+  return;
+end
+CPD = init_fields;
+ 
+CPD = class(CPD, 'gaussian_CPD', generic_CPD(0));
+
+args = varargin;
+ns = bnet.node_sizes;
+ps = parents(bnet.dag, self);
+dps = myintersect(ps, bnet.dnodes);
+cps = myintersect(ps, bnet.cnodes);
+fam_sz = ns([ps self]);
+
+CPD.self = self;
+CPD.sizes = fam_sz;
+
+% Figure out which (if any) of the parents are discrete, and which cts, and how big they are
+% dps = discrete parents, cps = cts parents
+CPD.cps = find_equiv_posns(cps, ps); % cts parent index
+CPD.dps = find_equiv_posns(dps, ps);
+ss = fam_sz(end);
+psz = fam_sz(1:end-1);
+dpsz = prod(psz(CPD.dps));
+cpsz = sum(psz(CPD.cps));
+
+% set default params
+CPD.mean = randn(ss, dpsz);
+CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]);    
+CPD.weights = randn(ss, cpsz, dpsz);
+CPD.cov_type = 'full';
+CPD.tied_cov = 0;
+CPD.clamped_mean = 0;
+CPD.clamped_cov = 0;
+CPD.clamped_weights = 0;
+CPD.cov_prior_weight = 0.01;
+CPD.cov_prior_entropic = 0;
+nargs = length(args);
+if nargs > 0
+  CPD = set_fields(CPD, args{:});
+end
+
+% Make sure the matrices have 1 dimension per discrete parent.
+% Bug fix due to Xuejing Sun 3/6/01
+CPD.mean = myreshape(CPD.mean, [ss ns(dps)]);
+CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]);
+CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]);
+
+% Precompute indices into block structured  matrices
+% to speed up CPD_to_lambda_msg and CPD_to_pi
+cpsizes = CPD.sizes(CPD.cps);
+CPD.cps_block_ndx = cell(1, length(cps));
+for i=1:length(cps)
+  CPD.cps_block_ndx{i} = block(i, cpsizes);
+end
+
+%%%%%%%%%%% 
+% Learning stuff
+
+% expected sufficient statistics 
+CPD.Wsum = zeros(dpsz,1);
+CPD.WYsum = zeros(ss, dpsz);
+CPD.WXsum = zeros(cpsz, dpsz);
+CPD.WYYsum = zeros(ss, ss, dpsz);
+CPD.WXXsum = zeros(cpsz, cpsz, dpsz);
+CPD.WXYsum = zeros(cpsz, ss, dpsz);
+
+% For BIC
+CPD.nsamples = 0;
+switch CPD.cov_type
+ case 'full',
+  % since symmetric 
+    %ncov_params = ss*(ss-1)/2; 
+    ncov_params = ss*(ss+1)/2; 
+  case 'diag',
+    ncov_params = ss;
+  otherwise
+    error(['unrecognized cov_type ' cov_type]);
+end
+% params = weights + mean + cov
+if CPD.tied_cov
+  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params;
+else
+  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params;
+end
+
+% for speeding up maximize_params
+CPD.useC = exist('rep_mult');
+
+clamped = CPD.clamped_mean && CPD.clamped_cov && CPD.clamped_weights;
+CPD = set_clamped(CPD, clamped);
+
+%%%%%%%%%%%
+
+function CPD = init_fields()
+% This ensures we define the fields in the same order 
+% no matter whether we load an object from a file,
+% or create it from scratch. (Matlab requires this.)
+
+CPD.self = [];
+CPD.sizes = [];
+CPD.cps = [];
+CPD.dps = [];
+CPD.mean = [];
+CPD.cov = [];
+CPD.weights = [];
+CPD.clamped_mean = [];
+CPD.clamped_cov = [];
+CPD.clamped_weights = [];
+CPD.cov_type = [];
+CPD.tied_cov = [];
+CPD.Wsum = [];
+CPD.WYsum = [];
+CPD.WXsum = [];
+CPD.WYYsum = [];
+CPD.WXXsum = [];
+CPD.WXYsum = [];
+CPD.nsamples = [];
+CPD.nparams = [];            
+CPD.cov_prior_weight = [];
+CPD.cov_prior_entropic = [];
+CPD.useC = [];
+CPD.cps_block_ndx = [];
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/gaussian_CPD_params_given_dps.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/gaussian_CPD_params_given_dps.m
new file mode 100644
index 00000000..72231a76
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/gaussian_CPD_params_given_dps.m
@@ -0,0 +1,28 @@
+function [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence)
+% GAUSSIAN_CPD_PARAMS_GIVEN_EV_ON_DPS Extract parameters given evidence on all discrete parents
+% function [m, C, W] = gaussian_CPD_params_given_ev_on_dps(CPD, domain, evidence)
+
+ps = domain(1:end-1);
+dps = ps(CPD.dps);
+if isempty(dps)
+  m = CPD.mean;
+  C = CPD.cov;
+  W = CPD.weights;
+else
+  odom = domain(~isemptycell(evidence(domain)));
+  dops = myintersect(dps, odom);
+  dpvals = cat(1, evidence{dops});
+  if length(dops) == length(dps)
+    dpsizes = CPD.sizes(CPD.dps);
+    dpval = subv2ind(dpsizes, dpvals(:)');
+    m = CPD.mean(:, dpval);
+    C = CPD.cov(:, :, dpval);
+    W = CPD.weights(:, :, dpval);
+  else
+    map = find_equiv_posns(dops, dps);
+    index = mk_multi_index(length(dps), map, dpvals);
+    m = CPD.mean(:, index{:});
+    C = CPD.cov(:, :, index{:});
+    W = CPD.weights(:, :, index{:});
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/get_field.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/get_field.m
new file mode 100644
index 00000000..2a50e1ac
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/get_field.m
@@ -0,0 +1,19 @@
+function val = get_params(CPD, name)
+% GET_PARAMS Get the parameters (fields) for a gaussian_CPD object
+% val = get_params(CPD, name)
+%
+% The following fields can be accessed
+%
+% mean       - mu(:,i) is the mean given Q=i
+% cov        - Sigma(:,:,i) is the covariance given Q=i 
+% weights    - W(:,:,i) is the regression matrix given Q=i 
+%
+% e.g., mean = get_params(CPD, 'mean')
+
+switch name
+ case 'mean',      val = CPD.mean;
+ case 'cov',       val = CPD.cov;
+ case 'weights',   val = CPD.weights;
+ otherwise,
+  error(['invalid argument name ' name]);
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/learn_params.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/learn_params.m
new file mode 100644
index 00000000..7ae5cb52
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/learn_params.m
@@ -0,0 +1,31 @@
+function CPD = learn_params(CPD, fam, data, ns, cnodes)
+%function CPD = learn_params(CPD, fam, data, ns, cnodes)
+% LEARN_PARAMS Compute the maximum likelihood estimate of the params of a gaussian CPD given complete data
+% CPD = learn_params(CPD, fam, data, ns, cnodes)
+%
+% data(i,m) is the value of node i in case m (can be cell array).
+% We assume this node has a maximize_params method.
+
+ncases = size(data, 2);
+CPD = reset_ess(CPD);
+% make a fully observed joint distribution over the family
+fmarginal.domain = fam;
+fmarginal.T = 1;
+fmarginal.mu = [];
+fmarginal.Sigma = [];
+if ~iscell(data)
+  cases = num2cell(data);
+else
+  cases = data;
+end
+hidden_bitv = zeros(1, max(fam));
+for m=1:ncases
+  % specify (as a bit vector) which elements in the family domain are hidden
+  hidden_bitv = zeros(1, max(fmarginal.domain));
+  ev = cases(:,m);
+  hidden_bitv(find(isempty(ev)))=1;
+  CPD = update_ess(CPD, fmarginal, ev, ns, cnodes, hidden_bitv);
+end
+CPD = maximize_params(CPD);
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/log_prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/log_prob_node.m
new file mode 100644
index 00000000..ac10f8a3
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/log_prob_node.m
@@ -0,0 +1,49 @@
+function L = log_prob_node(CPD, self_ev, pev)
+% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (gaussian)
+% L = log_prob_node(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+
+if iscell(self_ev), usecell = 1; else usecell = 0; end
+
+use_log = 1;
+ncases = length(self_ev);
+nparents = length(CPD.sizes)-1;
+assert(ncases == size(pev, 2));
+
+if ncases == 0
+  L = 0;
+  return;
+end
+
+L = 0;
+for m=1:ncases
+  if isempty(CPD.dps)
+    i = 1;
+  else
+    if usecell
+      dpvals = cat(1, pev{CPD.dps, m});
+    else
+      dpvals = pev(CPD.dps, m);
+    end
+    i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)');
+  end
+  if usecell
+    y = self_ev{m};
+  else
+    y = self_ev(m);
+  end
+  if length(CPD.cps) == 0 
+    L = L + gaussian_prob(y, CPD.mean(:,i), CPD.cov(:,:,i), use_log);
+  else
+    if usecell
+      x = cat(1, pev{CPD.cps, m});
+    else
+      x = pev(CPD.cps, m);
+    end
+    L = L + gaussian_prob(y, CPD.mean(:,i) + CPD.weights(:,:,i)*x, CPD.cov(:,:,i), use_log);
+  end
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/maximize_params.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/maximize_params.m
new file mode 100644
index 00000000..1624cbf2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/maximize_params.m
@@ -0,0 +1,68 @@
+function CPD = maximize_params(CPD, temp)
+% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
+% CPD = maximize_params(CPD, temperature)
+%
+% Temperature is currently ignored.
+
+if ~adjustable_CPD(CPD), return; end
+
+
+if CPD.clamped_mean
+  cl_mean = CPD.mean;
+else
+  cl_mean = [];
+end
+
+if CPD.clamped_cov
+  cl_cov = CPD.cov;
+else
+  cl_cov = [];
+end
+
+if CPD.clamped_weights
+  cl_weights = CPD.weights;
+else
+  cl_weights = [];
+end
+
+[ssz psz Q] = size(CPD.weights);
+
+[ss cpsz dpsz] = size(CPD.weights); % ss = self size = ssz
+if cpsz > CPD.nsamples
+  fprintf('gaussian_CPD/maximize_params: warning: input dimension (%d) > nsamples (%d)\n', ...
+	  cpsz, CPD.nsamples);
+end
+
+prior =  repmat(CPD.cov_prior_weight*eye(ssz,ssz), [1 1 Q]);
+
+
+[CPD.mean, CPD.cov, CPD.weights] = ...
+    clg_Mstep(CPD.Wsum, CPD.WYsum, CPD.WYYsum, [], CPD.WXsum, CPD.WXXsum, CPD.WXYsum, ...
+	      'cov_type', CPD.cov_type, 'clamped_mean', cl_mean, ...
+	      'clamped_cov', cl_cov, 'clamped_weights', cl_weights, ...
+	      'tied_cov', CPD.tied_cov, ...
+	      'cov_prior', prior);
+
+if 0
+CPD.mean = reshape(CPD.mean, [ss dpsz]);
+CPD.cov = reshape(CPD.cov, [ss ss dpsz]);
+CPD.weights = reshape(CPD.weights, [ss cpsz dpsz]);
+end
+
+% Bug fix 11 May 2003 KPM
+% clg_Mstep collapses all discrete parents into one mega-node
+% but convert_to_CPT needs access to each parent separately
+sz = CPD.sizes;
+ss = sz(end);
+
+% Bug fix KPM 20 May 2003: 
+cpsz = sum(sz(CPD.cps));
+%if isempty(CPD.cps)
+%  cpsz = 0;
+%else
+%  cpsz = sz(CPD.cps);
+%end
+dpsz = sz(CPD.dps);
+CPD.mean = myreshape(CPD.mean, [ss dpsz]);
+CPD.cov = myreshape(CPD.cov, [ss ss dpsz]);
+CPD.weights = myreshape(CPD.weights, [ss cpsz dpsz]);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/maximize_params_debug.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/maximize_params_debug.m
new file mode 100644
index 00000000..a588756d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/maximize_params_debug.m
@@ -0,0 +1,189 @@
+function CPD = maximize_params(CPD, temp)
+% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
+% CPD = maximize_params(CPD, temperature)
+%
+% Temperature is currently ignored.
+
+if ~adjustable_CPD(CPD), return; end
+
+CPD1 = struct(new_maximize_params(CPD));
+CPD2 = struct(old_maximize_params(CPD));
+assert(approxeq(CPD1.mean, CPD2.mean))
+assert(approxeq(CPD1.cov, CPD2.cov))
+assert(approxeq(CPD1.weights, CPD2.weights))
+
+CPD = new_maximize_params(CPD);
+
+%%%%%%%
+function CPD = new_maximize_params(CPD)
+
+if CPD.clamped_mean
+  cl_mean = CPD.mean;
+else
+  cl_mean = [];
+end
+
+if CPD.clamped_cov
+  cl_cov = CPD.cov;
+else
+  cl_cov = [];
+end
+
+if CPD.clamped_weights
+  cl_weights = CPD.weights;
+else
+  cl_weights = [];
+end
+
+[ssz psz Q] = size(CPD.weights);
+
+prior =  repmat(CPD.cov_prior_weight*eye(ssz,ssz), [1 1 Q]);
+[CPD.mean, CPD.cov, CPD.weights] = ...
+    Mstep_clg('w', CPD.Wsum, 'YY', CPD.WYYsum, 'Y', CPD.WYsum, 'YTY', [], ...
+	      'XX', CPD.WXXsum, 'XY', CPD.WXYsum, 'X', CPD.WXsum, ...
+	      'cov_type', CPD.cov_type, 'clamped_mean', cl_mean, ...
+	      'clamped_cov', cl_cov, 'clamped_weights', cl_weights, ...
+	      'tied_cov', CPD.tied_cov, ...
+	      'cov_prior', prior);
+
+
+%%%%%%%%%%%
+
+function CPD = old_maximize_params(CPD)
+
+
+if ~adjustable_CPD(CPD), return; end
+
+%assert(approxeq(CPD.nsamples, sum(CPD.Wsum)));
+assert(~any(isnan(CPD.WXXsum)))
+assert(~any(isnan(CPD.WXYsum)))
+assert(~any(isnan(CPD.WYYsum)))
+
+[self_size cpsize dpsize] = size(CPD.weights);
+
+% Append 1s to the parents, and derive the corresponding cross products.
+% This is used when estimate the means and weights simultaneosuly,
+% and when estimatting Sigma.
+% Let x2 = [x 1]'
+XY = zeros(cpsize+1, self_size, 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(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' 
+for i=1:dpsize
+  XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y
+	       CPD.WYsum(:,i)']; % 1*Y
+  % [x  * [x' 1]  = [xx' x
+  %  1]              x'  1]
+  XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i);
+	       CPD.WXsum(:,i)'   CPD.Wsum(i)];
+  YY(:,:,i) = CPD.WYYsum(:,:,i);
+end
+
+w = CPD.Wsum(:);
+% Set any zeros to one before dividing
+% This is valid because w(i)=0 => WYsum(:,i)=0, etc
+w = w + (w==0);
+
+if CPD.clamped_mean
+  % Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent
+  % to estimating B and then adding the clamped_mean to the last column.
+  if ~CPD.clamped_weights
+    B = zeros(self_size, cpsize, dpsize);
+    for i=1:dpsize
+      if det(CPD.WXXsum(:,:,i))==0
+	B(:,:,i) = 0;
+      else
+	% Eqn 9 in table 2 of TR
+	%B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i));
+	B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))';
+      end
+    end
+    %CPD.weights = reshape(B, [self_size cpsize dpsize]);
+    CPD.weights = B;
+  end
+elseif CPD.clamped_weights % KPM 1/25/02
+  if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals
+    for i=1:dpsize
+      CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i);
+    end
+  end
+else % nothing is clamped, so estimate mean and weights simultaneously
+  B2 = zeros(self_size, cpsize+1, dpsize);
+  for i=1:dpsize
+    if det(XX(:,:,i))==0  % fix by U. Sondhauss 6/27/99
+      B2(:,:,i)=0;          
+    else                    
+      % Eqn 9 in table 2 of TR
+      %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
+      B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))';
+    end                   
+    CPD.mean(:,i) = B2(:,cpsize+1,i);
+    CPD.weights(:,:,i) = B2(:,1:cpsize,i);
+  end
+end
+
+% Let B2 = [W mu]
+if cpsize>0
+  B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]);
+end
+B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]);
+
+% To avoid singular covariance matrices,
+% we use the regularization method suggested in "A Quasi-Bayesian approach to estimating
+% parameters for mixtures of normal distributions", Hamilton 91.
+% If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma),
+% where gamma >=0 is a smoothing parameter (equivalent sample size of I prior)
+
+gamma = CPD.cov_prior_weight;
+
+if ~CPD.clamped_cov
+  if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99
+    Z = 1-temp;
+    % When temp > 1, Z is negative, so we are dividing by a smaller
+    % number, ie. increasing the variance.
+  else
+    Z = 0;
+  end
+  if CPD.tied_cov
+    S = zeros(self_size, self_size);
+    % Eqn 2 from table 2 in TR
+    for i=1:dpsize
+      S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
+    end
+    %denom = CPD.nsamples + gamma + Z;
+    denom = CPD.nsamples +  Z;
+    S = (S + gamma*eye(self_size)) / denom;
+    if strcmp(CPD.cov_type, 'diag')
+      S = diag(diag(S));
+    end
+    CPD.cov = repmat(S, [1 1 dpsize]);
+  else 
+    for i=1:dpsize      
+      % Eqn 1 from table 2 in TR
+      S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i);
+      %denom = w(i) + gamma + Z;
+      denom = w(i) + Z;
+      S = (S + gamma*eye(self_size)) / denom;
+      CPD.cov(:,:,i) = S;
+    end
+    if strcmp(CPD.cov_type, 'diag')
+      for i=1:dpsize      
+	CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i)));
+      end
+    end
+  end
+end
+
+
+check_covars = 0;
+min_covar = 1e-5;
+if check_covars % prevent collapsing to a point
+  for i=1:dpsize
+    if min(svd(CPD.cov(:,:,i))) < min_covar
+      disp(['resetting singular covariance for node ' num2str(CPD.self)]);
+      CPD.cov(:,:,i) = CPD.init_cov(:,:,i);
+    end
+  end
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CPD_to_linear_gaussian.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CPD_to_linear_gaussian.m
new file mode 100644
index 00000000..dfc0cccc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CPD_to_linear_gaussian.m
@@ -0,0 +1,19 @@
+function [mu, Sigma, W] = CPD_to_linear_gaussian(CPD, domain, ns, cnodes, evidence)
+
+ps = domain(1:end-1);
+dnodes = mysetdiff(1:length(ns), cnodes);
+dps = myintersect(ps, dnodes); % discrete parents
+
+if isempty(dps)
+  Q = 1;
+else
+  assert(~any(isemptycell(evidence(dps))));
+  dpvals = cat(1, evidence{dps});
+  Q = subv2ind(ns(dps), dpvals(:)');
+end
+
+mu = CPD.mean(:,Q);
+Sigma = CPD.cov(:,:,Q);
+W = CPD.weights(:,:,Q);
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Entries
new file mode 100644
index 00000000..afb40930
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Entries
@@ -0,0 +1,2 @@
+/CPD_to_linear_gaussian.m/1.1.1.1/Wed May 29 15:59:52 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Repository
new file mode 100644
index 00000000..8aa921a1
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@gaussian_CPD/private
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/private/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/reset_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/reset_ess.m
new file mode 100644
index 00000000..d27105f0
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/reset_ess.m
@@ -0,0 +1,11 @@
+function CPD = reset_ess(CPD)
+% RESET_ESS Reset the Expected Sufficient Statistics for a Gaussian CPD.
+% CPD = reset_ess(CPD)
+
+CPD.nsamples = 0;    
+CPD.Wsum = zeros(size(CPD.Wsum));
+CPD.WYsum = zeros(size(CPD.WYsum));
+CPD.WYYsum = zeros(size(CPD.WYYsum));
+CPD.WXsum = zeros(size(CPD.WXsum));
+CPD.WXXsum = zeros(size(CPD.WXXsum));
+CPD.WXYsum = zeros(size(CPD.WXYsum));
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/sample_node.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/sample_node.m
new file mode 100644
index 00000000..74875eeb
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/sample_node.m
@@ -0,0 +1,22 @@
+function y = sample_node(CPD, pev)
+% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (gaussian)
+% y = sample_node(CPD, parent_evidence)
+%
+% pev{i} is the value of the i'th parent (if there are any parents)
+% y is the sampled value (a scalar or vector)
+
+if length(CPD.dps)==0
+  i = 1;
+else
+  dpvals = cat(1, pev{CPD.dps});
+  i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)');
+end
+
+if length(CPD.cps) == 0 
+  y = gsamp(CPD.mean(:,i), CPD.cov(:,:,i), 1);
+else
+  pev = pev(:);
+  x = cat(1, pev{CPD.cps});
+  y = gsamp(CPD.mean(:,i) + CPD.weights(:,:,i)*x(:), CPD.cov(:,:,i), 1);
+end
+y = y(:);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/set_fields.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/set_fields.m
new file mode 100644
index 00000000..4c1aef22
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/set_fields.m
@@ -0,0 +1,43 @@
+function CPD = set_fields(CPD, varargin)
+% SET_PARAMS Set the parameters (fields) for a gaussian_CPD object
+% CPD = set_params(CPD, name/value pairs)
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+%
+% mean       - mu(:,i) is the mean given Q=i
+% cov        - Sigma(:,:,i) is the covariance given Q=i 
+% weights    - W(:,:,i) is the regression matrix given Q=i 
+% cov_type   - if 'diag', Sigma(:,:,i) is diagonal 
+% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i
+% clamp_mean - if 1, we do not adjust mu(:,i) during learning 
+% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning 
+% clamp_weights - if 1, we do not adjust W(:,:,i) during learning
+% clamp      - if 1, we do not adjust any params
+% cov_prior_weight - weight given to I prior for estimating Sigma
+% cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0]
+%
+% e.g., CPD = set_params(CPD, 'mean', [0;0])
+
+args = varargin;
+nargs = length(args);
+for i=1:2:nargs
+  switch args{i},
+   case 'mean',        CPD.mean = args{i+1}; 
+   case 'cov',         CPD.cov = args{i+1}; 
+   case 'weights',     CPD.weights = args{i+1}; 
+   case 'cov_type',    CPD.cov_type = args{i+1}; 
+   %case 'tied_cov',    CPD.tied_cov = strcmp(args{i+1}, 'yes');
+   case 'tied_cov',    CPD.tied_cov = args{i+1};
+   case 'clamp_mean',  CPD.clamped_mean = args{i+1};
+   case 'clamp_cov',   CPD.clamped_cov = args{i+1};
+   case 'clamp_weights',  CPD.clamped_weights = args{i+1};
+   case 'clamp',  clamp = args{i+1};
+    CPD.clamped_mean = clamp;
+    CPD.clamped_cov = clamp;
+    CPD.clamped_weights = clamp;
+   case 'cov_prior_weight',  CPD.cov_prior_weight = args{i+1};
+   case 'cov_prior_entropic',  CPD.cov_prior_entropic = args{i+1};
+   otherwise,  
+    error(['invalid argument name ' args{i}]);
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/update_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/update_ess.m
new file mode 100644
index 00000000..3b58c02e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/update_ess.m
@@ -0,0 +1,88 @@
+function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+% UPDATE_ESS Update the Expected Sufficient Statistics of a Gaussian node
+% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+
+%if nargin < 6
+%  hidden_bitv = zeros(1, max(fmarginal.domain));
+%  hidden_bitv(find(isempty(evidence)))=1;
+%end
+
+dom = fmarginal.domain;
+self = dom(end);
+ps = dom(1:end-1);
+cps = myintersect(ps, cnodes);
+dps = mysetdiff(ps, cps);
+
+CPD.nsamples = CPD.nsamples + 1;            
+[ss cpsz dpsz] = size(CPD.weights); % ss = self size
+[ss dpsz] = size(CPD.mean);
+
+% Let X be the cts parent (if any), Y be the cts child (self).
+
+if ~hidden_bitv(self) && ~any(hidden_bitv(cps)) && all(hidden_bitv(dps))
+  % Speedup for the common case that all cts nodes are observed, all discrete nodes are hidden
+  % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0
+  % Since discrete parents are hidden, we do not need to add evidence to w.
+  w = fmarginal.T(:);
+  CPD.Wsum = CPD.Wsum + w;
+  y = evidence{self};
+  Cyy = y*y';
+  if ~CPD.useC
+     WY = repmat(w(:)',ss,1); % WY(y,i) = w(i)
+     WYY = repmat(reshape(WY, [ss 1 dpsz]), [1 ss 1]); % WYY(y,y',i) = w(i)
+     %CPD.WYsum = CPD.WYsum +  WY .* repmat(y(:), 1, dpsz);
+     CPD.WYsum = CPD.WYsum +  y(:) * w(:)';
+     CPD.WYYsum = CPD.WYYsum + WYY  .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]);
+  else
+     W = w(:)';
+     W2 = reshape(W, [1 1 dpsz]);
+     CPD.WYsum = CPD.WYsum +  rep_mult(W, y(:), size(CPD.WYsum)); 
+     CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum));
+  end
+  if cpsz > 0 % X exists
+    x = cat(1, evidence{cps}); x = x(:);
+    Cxx = x*x';
+    Cxy = x*y';
+    WX = repmat(w(:)',cpsz,1); % WX(x,i) = w(i)
+    WXX = repmat(reshape(WX, [cpsz 1 dpsz]), [1 cpsz 1]); % WXX(x,x',i) = w(i)
+    WXY = repmat(reshape(WX, [cpsz 1 dpsz]), [1 ss 1]); % WXY(x,y,i) = w(i)
+    if ~CPD.useC
+      CPD.WXsum = CPD.WXsum + WX .* repmat(x(:), 1, dpsz);
+      CPD.WXXsum = CPD.WXXsum + WXX .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]);
+      CPD.WXYsum = CPD.WXYsum + WXY .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]);
+    else
+      CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum));
+      CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum));
+      CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum));
+    end
+  end
+  return;
+end
+
+% general (non-vectorized) case
+fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow!
+
+if dpsz == 1 % no discrete parents
+  w = 1;
+else
+  w = fullm.T(:);
+end
+
+CPD.Wsum = CPD.Wsum + w;
+xi = 1:cpsz;
+yi = (cpsz+1):(cpsz+ss);
+for i=1:dpsz
+  muY = fullm.mu(yi, i);
+  SYY = fullm.Sigma(yi, yi, i);
+  CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
+  CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
+  if cpsz > 0
+    muX = fullm.mu(xi, i);
+    SXX = fullm.Sigma(xi, xi, i);
+    SXY = fullm.Sigma(xi, yi, i);
+    CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
+    CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
+    CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
+  end
+end                
+