Merge pull request #1 from ziejd2/octave_php_separate

Octave php separate

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(:);