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

14 files changed, 332 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CPD_to_lambda_msg.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CPD_to_lambda_msg.m
new file mode 100644
index 00000000..c323e8e5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CPD_to_lambda_msg.m
@@ -0,0 +1,62 @@
+function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
+% CPD_TO_LAMBDA_MSG Compute lambda message (gmux)
+% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
+% Pearl p183 eq 4.52
+
+% Let Y be this node, X1..Xn be the cts parents and M the discrete switch node.
+% e.g., for n=3, M=1
+%
+%  X1 X2 X3 M
+%   \
+%    \
+%      Y
+%
+% So the only case in which we send an informative message is if p=1=M.
+% To the other cts parents, we send the "know nothing" message.
+
+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);
+  dps = ps(CPD.dps);
+  M = evidence{dps};
+  if isempty(M)
+    error('gmux node must have observed discrete parent')
+  end
+  P = msg{n}.lambda.precision;
+  if all(P == 0) | (cps(M) ~= p) % if we know nothing, or are sending to a disconnected parent
+    lam_msg.precision = zeros(psz, psz);
+    lam_msg.info_state = zeros(psz, 1);
+    return;
+  end
+  % We are sending a message to the only effectively connected parent.
+  % There are no other incoming pi messages.
+  Bmu = CPD.mean(:,M);
+  BSigma = CPD.cov(:,:,M);
+  Bi = CPD.weights(:,:,M);
+  if (det(P) > 0) | 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
+    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/@gmux_CPD/CPD_to_pi.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CPD_to_pi.m
new file mode 100644
index 00000000..63b5726b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CPD_to_pi.m
@@ -0,0 +1,18 @@
+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',
+  dps = ps(CPD.dps);
+  k = evidence{dps};
+  if isempty(k)
+    error('gmux node must have observed discrete parent')
+  end
+  m = msg{n}.pi_from_parent{k}; 
+  B = CPD.weights(:,:,k);
+  pi.mu = CPD.mean(:,k) + B * m.mu;
+  pi.Sigma = CPD.cov(:,:,k) + B * m.Sigma * B';
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Entries
new file mode 100644
index 00000000..2a911068
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Entries
@@ -0,0 +1,7 @@
+/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:54 2002//
+/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/display.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/gmux_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/sample_node.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Entries.Log
new file mode 100644
index 00000000..24f16336
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/Old////
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Repository
new file mode 100644
index 00000000..8d764710
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@gmux_CPD
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Entries
new file mode 100644
index 00000000..f5a137ab
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Entries
@@ -0,0 +1,2 @@
+/gmux_CPD.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Repository
new file mode 100644
index 00000000..20395ac5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@gmux_CPD/Old
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/gmux_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/gmux_CPD.m
new file mode 100644
index 00000000..5c9507cf
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/Old/gmux_CPD.m
@@ -0,0 +1,92 @@
+function CPD = gmux_CPD(bnet, self, varargin)
+% GMUX_CPD Make a Gaussian multiplexer node
+%
+% CPD = gmux_CPD(bnet, node, ...) is used similarly to gaussian_CPD,
+% except we assume there is exactly one discrete parent (call it M)
+% which is used to select which cts parent to pass through to the output.
+% i.e., we define P(Y=y|M=m, X1, ..., XK) = N(y | W*x(m) + mu, Sigma)
+% where Y represents this node, and the Xi's are the cts parents.
+% All the Xi must have the same size, and the num values for M must be K.
+%
+% Currently the params for this kind of CPD cannot be learned.
+%
+% Optional arguments [ default in brackets ]
+%
+% mean       - mu  [zeros(Y,1)]
+% cov        - Sigma [eye(Y,Y)]
+% weights    - W [ randn(Y,X) ]
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  clamp = 0;
+  CPD = class(CPD, 'gmux_CPD', generic_CPD(clamp));
+  return;
+elseif isa(bnet, 'gmux_CPD')
+  % This might occur if we are copying an object.
+  CPD = bnet;
+  return;
+end
+CPD = init_fields;
+ 
+CPD = class(CPD, 'gmux_CPD', generic_CPD(1));
+
+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);
+if length(CPD.dps) ~= 1
+  error('gmux must have exactly 1 discrete parent')
+end
+ss = fam_sz(end);
+cpsz = fam_sz(CPD.cps(1)); % in gaussian_CPD, cpsz = sum(fam_sz(CPD.cps))
+if ~all(fam_sz(CPD.cps) == cpsz)
+  error('all cts parents must have same size')
+end
+dpsz = fam_sz(CPD.dps);
+if dpsz ~= length(cps)
+  error(['the arity of the mux node is ' num2str(dpsz) ...
+	 ' but there are ' num2str(length(cps)) ' cts parents']);
+end
+
+% set default params
+CPD.mean = zeros(ss, 1);
+CPD.cov = eye(ss);
+CPD.weights = randn(ss, cpsz);
+
+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}; 
+   otherwise,  
+    error(['invalid argument name ' args{i}]);
+  end
+end
+
+%%%%%%%%%%%
+
+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 = [];
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/convert_to_pot.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/convert_to_pot.m
new file mode 100644
index 00000000..bf8c29c4
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/convert_to_pot.m
@@ -0,0 +1,37 @@
+function pot = convert_to_pot(CPD, pot_type, domain, evidence)
+% CONVERT_TO_POT Convert a gmux CPD to a Gaussian potential
+% pot = convert_to_pot(CPD, pot_type, domain, evidence)
+  
+switch pot_type
+ case {'d', 'u', 'cg', 'scg'},
+  error(['can''t convert gmux to potential of type ' pot_type])
+
+ case {'c','g'},
+  % We create a large weight matrix with zeros in all blocks corresponding
+  % to the non-chosen parents, since they are effectively disconnected.
+  % The chosen parent is determined by the value, m,  of the discrete parent.
+  % Thus the potential is as large as the whole family.
+  ps = domain(1:end-1);
+  dps = ps(CPD.dps); % CPD.dps is an index, not a node number (because of param tying)
+  cps = ps(CPD.cps);
+  m = evidence{dps};
+  if isempty(m)
+    error('gmux node must have observed discrete parent')
+  end
+  bs = CPD.sizes(CPD.cps);
+  b = block(m, bs);
+  sum_cpsz = sum(CPD.sizes(CPD.cps));
+  selfsz = CPD.sizes(end);
+  W = zeros(selfsz, sum_cpsz);
+  W(:,b) = CPD.weights(:,:,m);
+
+  ns = zeros(1, max(domain));
+  ns(domain) = CPD.sizes;
+  self = domain(end);
+  cdom = [cps(:)' self];
+  pot = linear_gaussian_to_cpot(CPD.mean(:,m), CPD.cov(:,:,m), W, domain, ns, cdom, evidence);
+  
+ otherwise,
+  error(['unrecognized pot_type' pot_type])
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/display.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/display.m
new file mode 100644
index 00000000..4b04168c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/display.m
@@ -0,0 +1,4 @@
+function display(CPD)
+
+disp('gmux_CPD object');
+disp(struct(CPD));
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/gmux_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/gmux_CPD.m
new file mode 100644
index 00000000..4cef195c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/gmux_CPD.m
@@ -0,0 +1,95 @@
+function CPD = gmux_CPD(bnet, self, varargin)
+% GMUX_CPD Make a Gaussian multiplexer node
+%
+% CPD = gmux_CPD(bnet, node, ...) is used similarly to gaussian_CPD,
+% except we assume there is exactly one discrete parent (call it M)
+% which is used to select which cts parent to pass through to the output.
+% i.e., we define P(Y=y|M=m, X1, ..., XK) = N(y | W(m)*x(m) + mu(m), Sigma(m))
+% where Y represents this node, and the Xi's are the cts parents.
+% All the Xi must have the same size, and the num values for M must be K.
+%
+% Currently the params for this kind of CPD cannot be learned.
+%
+% Optional arguments [ default in brackets ]
+%
+% mean       - mu(:,i) is the mean given M=i [ zeros(Y,K) ]
+% cov        - Sigma(:,:,i) is the covariance given M=i [ repmat(1*eye(Y,Y), [1 1 K]) ]
+% weights    - W(:,:,i) is the regression matrix given M=i [ randn(Y,X,K) ]
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  clamp = 0;
+  CPD = class(CPD, 'gmux_CPD', generic_CPD(clamp));
+  return;
+elseif isa(bnet, 'gmux_CPD')
+  % This might occur if we are copying an object.
+  CPD = bnet;
+  return;
+end
+CPD = init_fields;
+ 
+CPD = class(CPD, 'gmux_CPD', generic_CPD(1));
+
+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);
+if length(CPD.dps) ~= 1
+  error('gmux must have exactly 1 discrete parent')
+end
+ss = fam_sz(end);
+cpsz = fam_sz(CPD.cps(1)); % in gaussian_CPD, cpsz = sum(fam_sz(CPD.cps))
+if ~all(fam_sz(CPD.cps) == cpsz)
+  error('all cts parents must have same size')
+end
+dpsz = fam_sz(CPD.dps);
+if dpsz ~= length(cps)
+  error(['the arity of the mux node is ' num2str(dpsz) ...
+	 ' but there are ' num2str(length(cps)) ' cts parents']);
+end
+
+% set default params
+%CPD.mean = zeros(ss, 1);
+%CPD.cov = eye(ss);
+%CPD.weights = randn(ss, cpsz);
+CPD.mean = zeros(ss, dpsz);
+CPD.cov = 1*repmat(eye(ss), [1 1 dpsz]);    
+CPD.weights = randn(ss, cpsz, dpsz);
+
+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}; 
+   otherwise,  
+    error(['invalid argument name ' args{i}]);
+  end
+end
+
+%%%%%%%%%%%
+
+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 = [];
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/sample_node.m b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/sample_node.m
new file mode 100644
index 00000000..53842a5d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gmux_CPD/sample_node.m
@@ -0,0 +1,10 @@
+function y = sample_node(CPD, pev)
+% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (gmux)
+% y = sample_node(CPD, parent_evidence)
+%
+% parent_ev{i} is the value of the i'th parent
+
+dpval = pev{CPD.dps};
+x = pev{CPD.cps(dpval)};
+y = gsamp(CPD.mean(:,dpval) + CPD.weights(:,:,dpval)*x(:), CPD.cov(:,:,dpval), 1);
+y = y(:);