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

18 files changed, 560 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries
new file mode 100644
index 00000000..1e0984dc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries
@@ -0,0 +1,11 @@
+/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/convert_to_table.m/1.1.1.1/Tue Mar 30 17:19:22 2004//
+/display.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/get_field.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/maximize_params.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/reset_ess.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//
+/set_fields.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/softmax_CPD.m/1.1.1.1/Tue Jan  7 16:25:14 2003//
+/update_ess.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries.Log
new file mode 100644
index 00000000..b2cd71e0
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Entries.Log
@@ -0,0 +1 @@
+A D/private////
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Repository
new file mode 100644
index 00000000..d5dac28b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@softmax_CPD
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_pot.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_pot.m
new file mode 100644
index 00000000..518f4a50
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_pot.m
@@ -0,0 +1,58 @@
+function pot = convert_to_pot(CPD, pot_type, domain, evidence)
+% CONVERT_TO_POT Convert a softmax CPD to a potential
+% pots = convert_to_pot(CPD, pot_type, domain, evidence)
+%
+% pots = CPD evaluated using evidence(domain)
+
+ncases = size(domain,2);
+assert(ncases==1); % not yet vectorized
+
+sz = dom_sizes(CPD);
+ns = zeros(1, max(domain));
+ns(domain) = sz;
+
+odom = domain(~isemptycell(evidence(domain)));
+T = convert_to_table(CPD, domain, evidence);
+
+switch pot_type
+ case 'u',
+  pot = upot(domain, sz, T, 0*myones(sz));  
+ case 'd',
+  ns(odom) = 1;
+  pot = dpot(domain, ns(domain), T);          
+ 
+ case {'c','g'},
+  % Since we want the output to be a Gaussian, the whole family must be observed.
+  % In other words, the potential is really just a constant.
+  p = T;
+  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));
+  ns(domain) = 0;
+  pot = cpot(domain, ns(domain), log(p));       
+ 
+ case 'cg',
+  T = T(:);
+  ns(odom) = 1;
+  can = cell(1, length(T));
+  for i=1:length(T)
+    can{i} = cpot([], [], log(T(i)));
+  end
+  ps = domain(1:end-1);
+  dps = ps(CPD.dpndx);
+  cps = ps(CPD.cpndx);
+  ddom = [dps CPD.self];
+  cdom = cps;
+  pot = cgpot(ddom, cdom, ns, can);   
+  
+ case 'scg'
+  T = T(:);
+  ns(odom) = 1;
+  pot_array = cell(1, length(T));
+  for i=1:length(T)
+    pot_array{i} = scgcpot([], [], T(i));
+  end
+  pot = scgpot(domain, [], [], ns, pot_array);   
+
+ otherwise,
+  error(['unrecognized pot type ' pot_type])
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_table.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_table.m
new file mode 100644
index 00000000..f703d79b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/convert_to_table.m
@@ -0,0 +1,52 @@
+function T = convert_to_table(CPD, domain, evidence)
+% CONVERT_TO_TABLE Convert a softmax CPD to a table, incorporating any evidence 
+% T = convert_to_table(CPD, domain, evidence)
+
+self       = domain(end);             
+ps         = domain(1:end-1);                            
+cnodes     = domain(CPD.cpndx);
+cps        = myintersect(ps, cnodes);
+dps        = domain(CPD.dpndx); 
+dps_as_cps = domain(CPD.dps_as_cps.ndx);
+all_dps    = union(dps,dps_as_cps);
+odom       = domain(~isemptycell(evidence(domain))); 
+if ~isempty(cps), assert(myismember(cps, odom)); end % all cts parents must be observed
+
+ns         = zeros(1, max(domain));
+ns(domain) = CPD.sizes;
+ens        = ns; % effective node sizes
+ens(odom)  = 1;
+
+% dpsize >= glimsz because the glm parameters are tied across the dps_as_cps parents
+dpsize       = prod(ens(all_dps)); % size of ALL self'discrete parents
+dpvals       = cat(1, evidence{myintersect(all_dps, odom)});
+cpvals       = cat(1, evidence{cps});
+if ~isempty(dps_as_cps),
+  separator          = CPD.dps_as_cps.separator;
+  dp_as_cpmap        = find_equiv_posns(dps_as_cps, all_dps);
+  dops_map           = find_equiv_posns(myintersect(all_dps, odom), all_dps);
+  puredp_map         = find_equiv_posns(dps, all_dps);
+  subs               = ind2subv(ens(all_dps), 1:prod(ens(all_dps)));
+  if ~isempty(dops_map), subs(:,dops_map) = subs(:,dops_map)+repmat(dpvals(:)',[size(subs,1) 1])-1; end
+end
+
+[w,b] = extract_params(CPD);
+T = zeros(dpsize, ns(self));                                       
+for i=1:dpsize,    
+  active_glm  = i;
+  dp_as_cpvals=zeros(1,sum(ns(dps_as_cps)));                                                                  
+  if ~isempty(dps_as_cps),                          
+    active_glm = max([1,subv2ind(ns(dps), subs(i,puredp_map))]);
+    % Extract the params compatible with the observations (if any) on the 'pure' discrete parents (if any)
+    where_one = separator + subs(i,dp_as_cpmap);
+    % and get in the dp_as_cp parents...
+    dp_as_cpvals(where_one)=1;                    
+  end                                               
+  T(i,:) = normalise(exp([dp_as_cpvals(:); cpvals(:)]'*w(:,:,active_glm) + b(:,active_glm)'));
+end
+if myismember(self, odom)
+  r = evidence{self};
+  T = T(:,r);
+end
+
+T = myreshape(T, ens(domain));
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/display.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/display.m
new file mode 100644
index 00000000..06a0f02c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/display.m
@@ -0,0 +1,4 @@
+function display(CPD)
+
+disp('softmax_CPD object');
+disp(struct(CPD)); 
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/get_field.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/get_field.m
new file mode 100644
index 00000000..240f1fd7
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/get_field.m
@@ -0,0 +1,18 @@
+function val = get_params(CPD, name)
+% GET_PARAMS Get the parameters (fields) for a softmax_CPD object
+% val = get_params(CPD, name)
+%
+% The following fields can be accessed
+%
+% weights - W(X,Y,Q)
+% offset  - b(Y,Q)
+%
+% e.g., W = get_params(CPD, 'weights')
+
+[W, b] = extract_params(CPD);
+switch name
+ case 'weights',   val = W;
+ case 'offset',    val = b;
+ otherwise,
+  error(['invalid argument name ' name]);
+end                
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/maximize_params.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/maximize_params.m
new file mode 100644
index 00000000..15c94dd5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/maximize_params.m
@@ -0,0 +1,41 @@
+function CPD = maximize_params(CPD, temp)
+% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (dsoftmax) using IRLS
+% CPD = maximize_params(CPD, temperature)
+% temperature parameter is ignored
+
+% Written by Pierpaolo Brutti
+
+if ~adjustable_CPD(CPD), return; end
+options = foptions;
+
+if CPD.verbose
+  options(1) = 1;
+else
+  options(1) = -1;
+end
+%options(1) = CPD.verbose;
+
+options(2) = CPD.wthresh;
+options(3) = CPD.llthresh;
+options(5) = CPD.approx_hess;
+options(14) = CPD.max_iter;
+
+dpsize = size(CPD.self_vals,3);
+for i=1:dpsize,
+  mask=find(CPD.eso_weights(:,:,i)>0); % for adapting the parameters we use only positive weighted example
+  if  ~isempty(mask),
+    if ~isempty(CPD.dps_as_cps.ndx),
+        puredp_map = find_equiv_posns(CPD.dpndx, union(CPD.dpndx, CPD.dps_as_cps.ndx)); % find the glm  structure
+        subs       = ind2subv(CPD.sizes(union(CPD.dpndx, CPD.dps_as_cps.ndx)),i);       % that corrisponds to the
+        active_glm = max([1,subv2ind(CPD.sizes(CPD.dpndx), subs(puredp_map))]);         % i-th 'fictitious' example
+        
+        CPD.glim{active_glm} = netopt_weighted(CPD.glim{active_glm}, options, CPD.parent_vals(mask',:,i),...
+            CPD.self_vals(mask',:,i), CPD.eso_weights(mask',:,i), 'scg');
+    else
+        alfa = 0.4; if CPD.solo, alfa = 1; end % learning step = 1 <=> self is all alone in the net
+        CPD.glim{i} = glmtrain_weighted(CPD.glim{i}, options, CPD.parent_vals(mask',:),...
+            CPD.self_vals(mask',:,i), CPD.eso_weights(mask',:,i), alfa);
+    end               
+  end
+  mask=[];
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Entries
new file mode 100644
index 00000000..b6610f0d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Entries
@@ -0,0 +1,2 @@
+/extract_params.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Repository
new file mode 100644
index 00000000..1667449e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@softmax_CPD/private
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/extract_params.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/extract_params.m
new file mode 100644
index 00000000..486af06e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/private/extract_params.m
@@ -0,0 +1,18 @@
+function [W, b] = extract_params(CPD)
+
+% W(X,Y,Q), b(Y,Q)  where Y = ns(self), X = ns(cps), Q = prod(ns(dps))
+
+glimsz = prod(CPD.sizes(CPD.dpndx));
+ss = CPD.sizes(end);
+cpsz       = sum(CPD.sizes(CPD.cpndx));
+dp_as_cpsz = sum(CPD.sizes(CPD.dps_as_cps.ndx));
+W = zeros(dp_as_cpsz + cpsz, ss, glimsz);
+b = zeros(ss, glimsz);
+
+for i=1:glimsz
+  W(:,:,i) = CPD.glim{i}.w1;
+  b(:,i) = CPD.glim{i}.b1(:);
+end
+
+W = myreshape(W, [dp_as_cpsz + cpsz ss CPD.sizes(CPD.dpndx)]);
+b = myreshape(b, [ss CPD.sizes(CPD.dpndx)]);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/reset_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/reset_ess.m
new file mode 100644
index 00000000..abf7d54e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/reset_ess.m
@@ -0,0 +1,8 @@
+function CPD = reset_ess(CPD)
+% RESET_ESS Reset the Expected Sufficient Statistics for a CPD (dsoftmax)
+% CPD = reset_ess(CPD)
+
+CPD.parent_vals = [];
+CPD.eso_weights=[];
+CPD.self_vals = [];
+CPD.nsamples = 0;  
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/sample_node.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/sample_node.m
new file mode 100644
index 00000000..1c519049
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/sample_node.m
@@ -0,0 +1,14 @@
+function y = sample_node(CPD, pvals)
+% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (discrete)
+% y = sample_node(CPD, parent_evidence)
+%
+% parent_evidence{i} is the value of the i'th parent
+
+n = length(pvals)+1;
+dom = 1:n;
+%evidence = cell(1,n);
+%evidence(1:n-1) = pvals(:)';
+evidence = pvals;
+evidence{end+1} = [];
+T = convert_to_table(CPD, dom, evidence);
+y = sample_discrete(T);
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/set_fields.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/set_fields.m
new file mode 100644
index 00000000..6c64b197
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/set_fields.m
@@ -0,0 +1,45 @@
+function CPD = set_params(CPD, varargin)
+% SET_PARAMS Set the parameters (fields) for a softmax_CPD object
+% CPD = set_params(CPD, name/value pairs)
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y), Q1=ns(dps(1)), Q2=ns(dps(2)), ...
+%   where dps are the discrete parents; if there are no discrete parents, we set Q1=1.)
+%
+% weights - (W(:,j,a,b,...) - W(:,j',a,b,...)) is ppn to dec. boundary
+%           between j,j' given Q1=a,Q2=b,... [ randn(X,Y,Q1,Q2,...) ]
+% offset  - (offset(j,a,b,...) - offset(j',a,b,...)) is the offset to dec. boundary
+%           between j,j' given Q1=a,Q2=b,... [ randn(Y,Q1,Q2,...) ]
+% clamped     - 'yes' means don't adjust params during learning ['no']
+% max_iter    - the maximum number of steps to take [10]
+% verbose     - 'yes' means print the LL at each step of IRLS ['no']
+% wthresh     - convergence threshold for weights [1e-2]
+% llthresh    - convergence threshold for log likelihood [1e-2]
+% approx_hess - 'yes' means approximate the Hessian for speed ['no']
+%
+% e.g., CPD = set_params(CPD,'offset', zeros(ns(i),1));
+
+args = varargin;
+nargs = length(args);
+glimsz = prod(CPD.sizes(CPD.dpndx));
+for i=1:2:nargs
+  switch args{i},
+   case 'discrete',     str='nothing to do';   
+   case 'clamped',      CPD = set_clamped(CPD, strcmp(args{i+1}, 'yes'));
+   case 'max_iter',     CPD.max_iter = args{i+1};
+   case 'verbose',      CPD.verbose = strcmp(args{i+1}, 'yes');
+   case 'max_iter',     CPD.max_iter = args{i+1};
+   case 'wthresh',      CPD.wthresh = args{i+1};
+   case 'llthresh',     CPD.llthresh = args{i+1};
+   case 'approx_hess',  CPD.approx_hess = strcmp(args{i+1}, 'yes');
+   case 'weights',      for q=1:glimsz, CPD.glim{q}.w1 = args{i+1}(:,:,q); end; 
+   case 'offset',
+    if glimsz == 1
+      CPD.glim{1}.b1 = args{i+1};
+    else
+      for q=1:glimsz, CPD.glim{q}.b1 = args{i+1}(:,q); end; 
+    end
+   otherwise,  
+    error(['invalid argument name ' args{i}]);       
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/softmax_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/softmax_CPD.m
new file mode 100644
index 00000000..3d2e5153
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/softmax_CPD.m
@@ -0,0 +1,187 @@
+function CPD = softmax_CPD(bnet, self, varargin)
+% SOFTMAX_CPD Make a softmax (multinomial logit) CPD
+%
+% To define this CPD precisely, let W be an (m x n) matrix with W(i,:) = {i-th row of B} 
+% => we can define the following vectorial function:
+%    
+%                                   softmax: R^n |--> R^m  
+%                  softmax(z,i-th)=exp(W(i,:)*z)/sum_k(exp(W(k,:)*z))      
+%
+% (this constructor augments z with a one at the beginning to introduce an offset term (=bias, intercept))                                   
+% Now call the continuous (cts) and always observed (obs) parents X,
+% the discrete parents (if any) Q, and this node Y then we use the discrete parent(s) just  to index
+% the parameter vectors (c.f., conditional Gaussian nodes); that is:
+%                 prob(Y=i | X=x, Q=j) = softmax(x,i-th|j)
+% where '|j' means that we are using the j-th (m x n) parameters matrix W(:,:,j).
+% If there are no discrete parents, this is a regular softmax node.
+% If Y is binary, this is a logistic (sigmoid) function.
+%
+% CPD = softmax_CPD(bnet, node_num, ...) will create a softmax CPD with random parameters,
+% where node is the number of a node in this equivalence class.
+%
+% The following optional arguments can be specified in the form of name/value pairs:
+% [default value in brackets]
+% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y), Q1=ns(dps(1)), Q2=ns(dps(2)), ...
+% where dps are the discrete parents; if there are no discrete parents, we set Q1=1.)
+%
+% discrete - the discrete parents that we want to treat like the cts ones [ [] ]. 
+%            This can be used to define sigmoid belief network - see below the reference.             
+%            For example suppose that Y has one cts parents X and two discrete ones: Q, C1 where:
+%            -> Q is binary (1/2) and used just to index the parameters of 'self'
+%            -> C1 is ternary (1/2/3) and treated as a cts node <=> its values appear into the linear 
+%               part of the softmax function
+%            then:
+%                     prob(Y|X=x, Q=q, C1=c1)= softmax(W(:,:,q)' * y)
+%            where y = [1 | delta(C1,1) delta(C1,2) delta(C1,3) | x(:)']' and delta(Y,a)=indicator(Y=a).
+% weights - (w(:,j,a,b,...) - w(:,j',a,b,...)) is ppn to dec. boundary
+%           between j,j' given Q1=a,Q2=b,... [ randn(X,Y,Q1,Q2,...) ]
+% offset  - (b(j,a,b,...) - b(j',a,b,...)) is the offset to dec. boundary
+%           between j,j' given Q1=a,Q2=b,... [ randn(Y,Q1,Q2,...) ]
+%
+% e.g., CPD = softmax_CPD(bnet, i, 'offset', zeros(ns(i),1));
+%
+% The following fields control the behavior of the M step, which uses 
+% a weighted version of the Iteratively Reweighted Least Squares (WIRLS) if dps_as_cps=[]; or
+% a weighted SCG otherwise, as implemented in Netlab, and modified by Pierpaolo Brutti.
+%
+% clamped     - 'yes' means don't adjust params during learning ['no']
+% max_iter    - the maximum number of steps to take [10]
+% verbose     - 'yes' means print the LL at each step of IRLS ['no']
+% wthresh     - convergence threshold for weights [1e-2]
+% llthresh    - convergence threshold for log likelihood [1e-2]
+% approx_hess - 'yes' means approximate the Hessian for speed ['no']
+%
+% For backwards compatibility with BNT2, you can also specify the parameters in the following order
+%   softmax_CPD(bnet, self, w, b, clamped, max_iter, verbose, wthresh, llthresh, approx_hess)
+%
+% REFERENCE
+% For details on the sigmoid belief nets, see:
+% - Neal (1992). Connectionist learning of belief networks, Artificial Intelligence, 56, 71-113.
+% - Saul, Jakkola, Jordan (1996). Mean field theory for sigmoid belief networks, Journal of Artificial Intelligence Reseach (4), pagg. 61-76.
+%
+% For details on the M step, see:
+% - K. Chen, L. Xu, H. Chi (1999). Improved learning algorithms for mixtures of experts in multiclass 
+%       classification. Neural Networks 12, pp. 1229-1252.
+% - M.I. Jordan, R.A. Jacobs (1994). Hierarchical Mixtures of Experts and the EM algorithm. 
+%       Neural Computation 6, pp. 181-214.
+% - S.R. Waterhouse, A.J. Robinson (1994). Classification Using Hierarchical Mixtures of Experts. In Proc. IEEE
+%       Workshop on Neural Network for Signal Processing IV, pp. 177-186
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  CPD = class(CPD, 'softmax_CPD', discrete_CPD(0, []));
+  return;
+elseif isa(bnet, 'softmax_CPD')
+  % This might occur if we are copying an object.
+  CPD = bnet;
+  return;
+end
+CPD = init_fields;
+
+assert(myismember(self, bnet.dnodes));
+ns = bnet.node_sizes;
+ps = parents(bnet.dag, self);
+dps = myintersect(ps, bnet.dnodes);
+cps = myintersect(ps, bnet.cnodes);
+
+clamped = 0;
+CPD = class(CPD, 'softmax_CPD', discrete_CPD(clamped, ns([ps self])));
+
+dps_as_cpssz = 0;
+dps_as_cps = [];
+% determine if any discrete parents are to be treated as cts
+if nargin >= 3 && isstr(varargin{1}) % might have passed in 'discrete'
+  for i=1:2:length(varargin)
+    if strcmp(varargin{i}, 'discrete')
+      dps_as_cps = varargin{i+1};
+      assert(myismember(dps_as_cps, dps));
+      dps = mysetdiff(dps, dps_as_cps);         % put out the dps treated as cts
+      CPD.dps_as_cps.ndx = find_equiv_posns(dps_as_cps, ps);
+      CPD.dps_as_cps.separator = [0 cumsum(ns(dps_as_cps(1:end-1)))]; % concatenated dps_as_cps dims separators
+      dps_as_cpssz = sum(ns(dps_as_cps));
+      break;
+    end
+  end
+end
+assert(~isempty(union(cps, dps_as_cps)));   % It have to be at least a cts or a dps_as_cps parents
+self_size = ns(self); 
+cpsz = sum(ns(cps));  
+glimsz = prod(ns(dps));
+CPD.dpndx = find_equiv_posns(dps, ps);  % it contains only the indeces of the 'pure' dps
+CPD.cpndx = find_equiv_posns(cps, ps);
+
+CPD.self  = self;
+CPD.solo  = (length(ns)<=2);
+CPD.sizes = bnet.node_sizes([ps self]);
+
+% set default params
+CPD.max_iter = 10;
+CPD.verbose = 0;
+CPD.wthresh = 1e-2;
+CPD.llthresh = 1e-2;
+CPD.approx_hess = 0;
+CPD.glim = cell(1,glimsz);
+for i=1:glimsz
+  CPD.glim{i} = glm(dps_as_cpssz + cpsz, self_size, 'softmax');
+end
+
+if nargin >= 3
+  args = varargin;
+  nargs = length(args);
+  if ~isstr(args{1})
+    %   softmax_CPD(bnet, self, w, b, clamped, max_iter, verbose, wthresh, llthresh, approx_hess)
+    if nargs >= 1 && ~isempty(args{1}), CPD = set_fields(CPD, 'weights', args{1}); end
+    if nargs >= 2 && ~isempty(args{2}), CPD = set_fields(CPD, 'offset', args{2});  end
+    if nargs >= 3 && ~isempty(args{3}), CPD = set_clamped(CPD, args{3});           end
+    if nargs >= 4 && ~isempty(args{4}), CPD.max_iter    = args{4}; end
+    if nargs >= 5 && ~isempty(args{5}), CPD.verbose     = args{5}; end
+    if nargs >= 6 && ~isempty(args{6}), CPD.wthresh     = args{6}; end
+    if nargs >= 7 && ~isempty(args{7}), CPD.llthresh   = args{7}; end
+    if nargs >= 8 && ~isempty(args{8}), CPD.approx_hess = args{8}; end
+  else
+    CPD = set_fields(CPD, args{:});
+  end
+end
+
+% sufficient statistics 
+% Since dsoftmax is not in the exponential family, we must store all the raw data.
+CPD.parent_vals = [];         % X(l,:) = value of cts parents in l'th example
+CPD.self_vals = [];           % Y(l,:) = value of self in l'th example
+
+CPD.eso_weights=[];           % weights used by the WIRLS algorithm
+
+% For BIC
+CPD.nsamples = 0;   
+if ~adjustable_CPD(CPD),
+   CPD.nparams=0;
+else
+   [W, b] = extract_params(CPD);
+   CPD.nparams= prod(size(W)) + prod(size(b));
+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.glim = {};
+CPD.self = [];
+CPD.solo = [];
+CPD.max_iter = [];
+CPD.verbose = [];
+CPD.wthresh = [];
+CPD.llthresh = [];
+CPD.approx_hess = [];
+CPD.sizes = [];
+CPD.parent_vals = [];
+CPD.eso_weights=[];
+CPD.self_vals = [];
+CPD.nsamples = [];
+CPD.nparams = [];
+CPD.dpndx = [];
+CPD.cpndx = [];
+CPD.dps_as_cps.ndx = [];
+CPD.dps_as_cps.separator = [];
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/update_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/update_ess.m
new file mode 100644
index 00000000..143c567c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@softmax_CPD/update_ess.m
@@ -0,0 +1,97 @@
+function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+% UPDATE_ESS Update the Expected Sufficient Statistics of a softmax node
+% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+%
+% fmarginal = overall posterior distribution of self and its parents
+% fmarginal(i1,i2...,ik,s)=prob(Pa1=i1,...,Pak=ik, self=s| X)
+% 
+% => 1) prob(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak) with prob(Pa1,...,Pak)=sum{s,fmarginal}
+%       [self estimation -> CPD.self_vals]
+% 	  2) prob(Pa1,...,Pak) [WIRLS weights -> CPD.eso_weights]
+%
+% Hidden_bitv is ignored
+
+% Written by Pierpaolo Brutti
+
+if ~adjustable_CPD(CPD), return; end
+
+domain     = fmarginal.domain;                              
+self       = domain(end);          
+ps         = domain(1:end-1);                                     
+cnodes     = domain(CPD.cpndx);
+cps        = myintersect(domain, cnodes);                     
+dps        = mysetdiff(ps, cps);                            
+dn_use     = dps;
+if isempty(evidence{self}) dn_use = [dn_use self]; end % if self is hidden we must consider its dimension  
+dps_as_cps = domain(CPD.dps_as_cps.ndx);
+odom       = domain(~isemptycell(evidence(domain))); 
+
+ns = zeros(1, max(domain));
+ns(domain) = CPD.sizes;     % CPD.sizes = bnet.node_sizes([ps self]);
+ens = ns;                   % effective node sizes
+ens(odom) = 1;              
+dpsize = prod(ns(dps));
+
+% Extract the params compatible with the observations (if any) on the discrete parents (if any)
+dops = myintersect(dps, odom);
+dpvals = cat(1, evidence{dops});
+
+subs = ind2subv(ens(dn_use), 1:prod(ens(dn_use)));
+dpmap = find_equiv_posns(dops, dn_use);
+if ~isempty(dpmap), subs(:,dpmap) = subs(:,dpmap)+repmat(dpvals(:)',[size(subs,1) 1])-1; end
+supportedQs = subv2ind(ns(dn_use), subs); subs=subs(1:prod(ens(dps)),1:length(dps));
+Qarity = prod(ns(dn_use));
+if isempty(dn_use), Qarity = 1; end   
+
+fullm.T              = zeros(Qarity, 1);
+fullm.T(supportedQs) = fmarginal.T(:);
+rs_dim = CPD.sizes;    rs_dim(CPD.cpndx) = 1;           %
+if ~isempty(evidence{self}), rs_dim(end)=1; end         % reshaping the marginal
+fullm.T              = reshape(fullm.T, rs_dim);        %
+
+% --------------------------------------------------------------------------------UPDATE--
+
+CPD.nsamples = CPD.nsamples + 1;
+
+% 1) observations vector -> CPD.parents_vals ---------------------------------------------
+cpvals = cat(1, evidence{cps});
+
+if ~isempty(dps_as_cps),   % ...get in the dp_as_cp parents... 
+    separator          = CPD.dps_as_cps.separator;
+    dp_as_cpmap        = find_equiv_posns(dps_as_cps, dps);       
+    for i=1:dpsize,
+        dp_as_cpvals=zeros(1,sum(ns(dps_as_cps)));
+        possible_vals = ind2subv(ns(dps),i);
+        ll=find(ismember(subs(:,dp_as_cpmap), possible_vals(dp_as_cpmap), 'rows')==1);   
+        if ~isempty(ll),
+            where_one = separator + possible_vals(dp_as_cpmap);
+            dp_as_cpvals(where_one)=1;                            
+        end
+        CPD.parent_vals(CPD.nsamples,:,i) = [dp_as_cpvals(:); cpvals(:)]';
+    end
+else
+    CPD.parent_vals(CPD.nsamples,:) = cpvals(:)';
+end
+
+% 2) weights vector -> CPD.eso_weights ----------------------------------------------------
+if isempty(evidence{self}),             % self is hidden
+    pesi=reshape(sum(fullm.T, length(rs_dim)),[dpsize,1]);
+else
+    pesi=reshape(fullm.T,[dpsize,1]);
+end
+assert(approxeq(sum(pesi),1));          % check
+
+% 3) estimate (if R is hidden) or recover (if R is obs) self'value-------------------------
+if isempty(evidence{self})                                  % P(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak)
+    r=reshape(mk_stochastic(fullm.T), [dpsize ns(self)]);   % matrix size: prod{j,ns(Paj)} x ns(self)      
+else
+    r = zeros(dpsize,ns(self));
+    for i=1:dpsize, if pesi(i)~=0, r(i,evidence{self}) = 1; end; end
+end
+for i=1:dpsize, if pesi(i)~=0, assert(approxeq(sum(r(i,:)),1)); end; end     % check
+
+% 4) save the previous values --------------------------------------------------------------
+for i=1:dpsize
+    CPD.eso_weights(CPD.nsamples,:,i)=pesi(i);
+    CPD.self_vals(CPD.nsamples,:,i) = r(i,:); 
+end