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

8 files changed, 922 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab.m b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab.m
new file mode 100644
index 00000000..5bc375fc
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab.m
@@ -0,0 +1,179 @@
+function [dag,best_score] =	learn_struct_gs_dtab(data, nodesizes, seeddag, varargin)
+%
+% LEARN_STRUCT_GS(data,seeddag)	learns a structure of Bayesian net by Greedy Search.
+% dag =	learn_struct_gs(data, nodesizes, seeddag)
+%
+% dag: the final structurre	matrix
+% Data : training data,	data(i,m) is the m obsevation of node i
+% Nodesizes: the size array	of different nodes
+% seeddag: given seed Dag for hill climbing, optional
+%
+% by Gang Li @ Deakin University (gli73@hotmail.com)
+%
+% -----------------------------------------------------
+%
+% Modified from	learn_struct_gs (SLP 1.3)
+% to learn structure of	BN with	tabular	nodes:
+%
+% 1) make use of score decomposibility 
+% 2) still to do
+%
+% by Darima	<darrimma@yahoo.com>, 27/12/2005
+%
+
+[N ncases] = size(data);
+if (nargin < 3 ) 
+	seeddag	= zeros(N,N); %	mk_rnd_dag(N); %call BNT function
+elseif ~acyclic(seeddag)
+	seeddag	= mk_rnd_dag(N); %zeros(N,N);
+end;
+
+% set default params (the same as in score_dags)
+for i=1:N
+  type{i} = 'tabular';
+  params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
+end
+scoring_fn = 'bayesian';
+discrete = 1:N;
+verbose	 = 'yes';
+
+% get params
+args = varargin;
+nargs =	length(args);
+if length(args)	> 0
+	if isstr(args{1})
+		for	i =	1:2:nargs
+			switch args{i}
+			case 'scoring_fn', scoring_fn =	args{i+1};
+            case 'type',       type = args{i+1}; 
+            case 'discrete',   discrete = args{i+1}; 
+            case 'params',     
+                if isempty(args{i+1}), params = cell(1,n); 
+                else params = args{i+1};  end;           
+			case 'verbose',	 verbose  =	strcmp(args{i+1},'yes');
+			end;
+		end;
+	end;
+end;
+
+done = 0;
+best_score = score_dags(data,nodesizes,	{seeddag},'scoring_fn',scoring_fn);
+
+it = 0;
+while ~done
+	[dags,op,nodes]	= mk_nbrs_of_dag(seeddag);
+	nbrs = length(dags);
+%------ DEBUG    
+    fprintf('DEBUG: nbrs = %6d\n',nbrs);
+    time_inloop = cputime;
+%------ DEND
+    
+%  	scores	= score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn);
+	scores = zeros(nbrs,1);
+	for	i =	1:nbrs
+		xj = nodes(i,2);	
+		ps_old = parents(seeddag, xj)';
+		ps_new = parents(dags{i}, xj)';
+        scor_old = score_family(xj, ps_old, type{xj}, scoring_fn, ...
+            nodesizes, discrete, data, params{xj});
+        scor_new = score_family(xj, ps_new, type{xj}, scoring_fn, ...
+            nodesizes, discrete, data, params{xj});
+		scores(i) =	best_score - scor_old +	scor_new;
+		if isequal(op{i},'rev')
+		    xi = nodes(i,1);
+			ps_old = parents(seeddag, xi);
+			ps_new = parents(dags{i}, xi);     
+            scor_old = score_family(xi, ps_old, type{xi}, scoring_fn, ...
+                nodesizes, discrete, data, params{xi});
+            scor_new = score_family(xi, ps_new, type{xi}, scoring_fn, ...
+                nodesizes, discrete, data, params{xi});
+			scores(i) =	scores(i) - scor_old +	scor_new;            
+		end
+    end
+    
+	max_score =	max(scores);
+	new	= find(scores == max_score );
+    
+%------ DEBUG
+    fprintf('       -> max_score = %7.5f\n',max_score);
+    fprintf('       -> find(scores == max_score): %d... of %d\n',...
+        new(1),length(new));
+%------ DEND
+	if ~isempty(new) & (max_score >	best_score)
+		p =	sample_discrete(normalise(ones(1, length(new))));
+		best_score = max_score;
+		seeddag	= dags{new(p)};
+	else
+		done = 1;
+	end;
+    
+    it = it+1;
+%------ DEBUG
+    time_inloop = cputime-time_inloop;
+    fprintf('       time = %12.5f\n',time_inloop);
+%------ DEND
+end;
+
+dag	= seeddag;
+
+outcount = 0; 
+best_score = score_dags(data,nodesizes,	{seeddag},'scoring_fn',scoring_fn);
+while outcount < 2
+	innercount = 0;
+	for	i=1:N
+		for	j=1:N
+		   if i==j,	continue;	 end;
+		   if seeddag(i,j) == 0	 % No edge i-->j, then try to add it
+			   tempdag = seeddag;
+			   tempdag(i,j)	= 1;
+			   if acyclic(tempdag)
+					temp_score = score_dags(data,nodesizes,	{tempdag},'scoring_fn',scoring_fn);
+					if temp_score >	best_score
+						seeddag	= tempdag;
+						best_score=	temp_score;
+						innercount = innercount	+1;
+					end;
+			   end
+		   else	 % exists edge i--j, then try reverse it or	remove it
+			   tempdag = seeddag;
+			   tempdag(i,j)	= 0; tempdag(j,i) =	1; 
+			   if acyclic(tempdag)
+				   temp_score =	score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+				   if temp_score > best_score
+					   seeddag = tempdag;
+					   best_score =	temp_score;
+					   innercount =	innercount +1;
+				   else
+					   tempdag = seeddag;
+					   tempdag(i,j)	= 0;
+					   temp_score =	score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+					   if temp_score > best_score
+						   seeddag = tempdag;
+						   best_score= temp_score;
+						   innercount =	innercount +1;
+					   end;
+				   end;
+			   else
+				   tempdag = seeddag;
+				   tempdag(i,j)=0;
+				   temp_score =	score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+				   if temp_score > best_score
+					   seeddag = tempdag;
+					   best_score= temp_score;
+					   innercount =	innercount +1;
+				   end;
+			   end;
+		   end;
+		end; % end for j
+	end; % end for i
+	if innercount == 0
+		outcount = outcount	+1;
+	end;
+end;  %	end	while
+
+%------ DEBUG
+fprintf('DEBUG: Number of iterations = %d\n',it);
+fprintf('DEBUG: Outcount = %d\n',outcount);
+%------ DEND
+dag	= seeddag;
+
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab_INFO.txt b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab_INFO.txt
new file mode 100644
index 00000000..dc076d78
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab_INFO.txt
@@ -0,0 +1,91 @@
+Objet:  	Re: [BayesNetToolbox] cache size in structure learning package (SLP)	
+De:  	"Darima Lamazhapova" <darrimma@yahoo.com>	
+Date:  	Mar 3 janvier 2006 14:13	
+to:  	BayesNetToolbox@yahoogroups.com	
+
+Hi, 
+I have not done extensive tests with different dataset and cache sizes.
+Mostly because I was interested in the case of small dataset,
+say dataset size 100. Even in this case seems that using cache 
+drasticly improves the performance of greedy search. Here are some
+numbers (alarm network, dataset 100, greedy search):
+    no cache     - ~11000 sec
+    cache  100  -  1395 sec
+    cache  200  -  1386 sec
+    cache  300  -  1428 sec
+    cache  500  -  1567 sec
+    cache 1000 -  1680 sec
+    cache 3000 -  2196 sec
+The reason why greedy search takes so much time without using
+cache, i think is that it does not make use of score 
+decomposibility.
+When I have modified the code so that it used decomposibility,
+ time required for the search decreased from 11000 to 1362 sec
+without using cache. Reimplemening score_family in C decreased this
+time to 552 sec; after reimplementing score_dag in C and
+removing outcount cycle from the learn_struct_gs code (I did not
+understand why it is necessary, can anyone explain what is
+the use of it, please?) the time required for GS dropped up to  69 sec.
+Implementation in C was done only for tabular nodes, Bayesian scoring
+function with default priors. If anyone interested I can submit the codes,
+although they are not very well tested.
+Darima
+
+Olivier Francois <olivier.francois@insa-rouen.fr> wrote:
+
+> Hello,
+> Can anyone give recommendations on choosing cache size?
+> When I used 500, calculations using learn_struct_gs2 for
+> ALARM network took about 20 min (dataset size was 100).
+> I thought it is a bit slow for 3 GHz computer with 1 Gb memory,
+> or am I wrong here?
+> I thought that increasing cache size might increase performance,
+> however the calculations took even longer.
+> I am kinda lost right now.
+> Any comments would be highly appreciated.
+> Darima
+>
+
+Hi,
+
+I have seen this phenomenon.
+In fact, for all the tests I have done, I advise you to use a cache of size between 200
+and 500.
+
+When the size is bigger the time spent to search if an entry already exists is quite
+similar to the time spend to recalculate the score, specialy if you have a small dataset
+(under 1500-2000 samples).
+Nevertheless, if you've got a huge dataset (5000 or more), it will be very advantageous
+to use the cache option.
+
+Moreover, I have seen that when you use a big cache (1000 or more), it do not speed up
+the computationnal time but you not spent a lot of extented time in using it.
+I think it is better to have a too big cache than a too small.
+
+If the cache is too small you often erase entries that will be recalculate later,
+specially if you have a lot of attributes, and the time spent will be equivalent or
+higher than if you have not used it.
+
+But suprisigly, even if you have a lot of nodes and a lot of samples (I have tested to
+25x10000 if I remember well), upgrading the size of the cache (1000 and more) do not
+seem to improve the computation time.
+
+- More tests are needed to be sure of that -
+
+I have not done a lot of tests, and this is only what I believed to see.
+Maybe it is better to take a cache of size 2500 or 5000 for a dataset of size 40*50000
+or 60*5000...?
+
+I also have remarked that, if you use a 'sparse matrix' for the cache inst ead of a
+stardart one, that causes a waste of time even if the cache matrix is, in fact, sparse.
+
+
+In your case, with 100 samples, you don't need too use the cache option.
+
+If you make (or have made for some others ?) more tests with different sample sizes and
+different numbers of attributes, I am interrested in getting some comments on your use
+of this function.
+
+
+Bonne fetes - Happy Hollidays - Felice Fiestas
+    Olivier F.
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabx.m b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabx.m
new file mode 100644
index 00000000..be624ee8
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabx.m
@@ -0,0 +1,181 @@
+function [dag,best_score] =	learn_struct_gs_dtabx(data, nodesizes, seeddag, varargin)
+%
+% LEARN_STRUCT_GS(data,seeddag)	learns a structure of Bayesian net by Greedy Search.
+% dag =	learn_struct_gs(data, nodesizes, seeddag)
+%
+% dag: the final structure matrix
+% Data : training data,	data(i,m) is the m obsevation of node i
+% Nodesizes: the size array	of different nodes
+% seeddag: given seed Dag for hill climbing, optional
+%
+% by Gang Li @ Deakin University (gli73@hotmail.com)
+%
+% -----------------------------------------------------
+%
+% Modified from	learn_struct_gs (SLP 1.3)
+% to learn structure of	BN with	tabular	nodes:
+%
+% 1) make use of score decomposibility 
+% 2) replace score_family with score_family_x.c that calculates 
+%    Bayesian score with default parameters (non-adjustable)
+%
+% by Darima	<darrimma@yahoo.com>, 28/12/2005
+%
+
+[N ncases] = size(data);
+if (nargin < 3 ) 
+	seeddag	= zeros(N,N); %	mk_rnd_dag(N); %call BNT function
+elseif ~acyclic(seeddag)
+	seeddag	= mk_rnd_dag(N); %zeros(N,N);
+end;
+
+% set default params (the same as in score_dags)
+for i=1:N
+  type{i} = 'tabular';
+  params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
+end
+scoring_fn = 'bayesian';
+discrete = 1:N;
+verbose	 = 'yes';
+
+% get params
+args = varargin;
+nargs =	length(args);
+if length(args)	> 0
+	if isstr(args{1})
+		for	i =	1:2:nargs
+			switch args{i}
+			case 'scoring_fn', scoring_fn =	args{i+1};
+            case 'type',       type = args{i+1}; 
+            case 'discrete',   discrete = args{i+1}; 
+            case 'params',     
+                if isempty(args{i+1}), params = cell(1,n); 
+                else params = args{i+1};  end;           
+			case 'verbose',	 verbose  =	strcmp(args{i+1},'yes');
+			end;
+		end;
+	end;
+end;
+
+%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+datax = data-ones(N,ncases);     %%+++++++++++++++++++++++ IMPORTANT!!!!
+%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+best_score = score_dags(data,nodesizes,	{seeddag},'scoring_fn',scoring_fn);
+done = 0;
+it = 0;
+
+while ~done
+	[dags,op,nodes]	= mk_nbrs_of_dag(seeddag);
+	nbrs = length(dags);
+%------ DEBUG    
+    fprintf('DEBUG: nbrs = %6d\n',nbrs);
+    time_inloop = cputime;
+%------ DEND
+%     scores	= score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn);
+	scores = zeros(nbrs,1);
+	for	i =	1:nbrs
+		xj = nodes(i,2);	
+		ps_old = parents(seeddag, xj)';
+		ps_new = parents(dags{i}, xj)';
+        scor_old = score_family_x([datax(ps_old,:);datax(xj,:)],...
+            [nodesizes(ps_old),nodesizes(xj)]);
+        scor_new = score_family_x([datax(ps_new,:);datax(xj,:)],...
+            [nodesizes(ps_new),nodesizes(xj)]);
+		scores(i) =	best_score - scor_old +	scor_new;
+		if isequal(op{i},'rev')
+		    xi = nodes(i,1);
+			ps_old = parents(seeddag, xi);
+			ps_new = parents(dags{i}, xi);
+            scor_old = score_family_x([datax(ps_old,:);datax(xi,:)],...
+                [nodesizes(ps_old),nodesizes(xi)]);
+            scor_new = score_family_x([datax(ps_new,:);datax(xi,:)],...
+                [nodesizes(ps_new),nodesizes(xi)]);
+			scores(i) =	scores(i) - scor_old +	scor_new;            
+		end
+    end    
+	max_score =	max(scores);
+	new	= find(scores == max_score );
+%------ DEBUG
+    fprintf('       -> max_score = %7.5f\n',max_score);
+    fprintf('       -> find(scores == max_score): %d... of %d\n',...
+        new(1),length(new));
+%------ DEND
+	if ~isempty(new) & (max_score >	best_score)
+		p =	sample_discrete(normalise(ones(1, length(new))));
+		best_score = max_score;
+		seeddag	= dags{new(p)};
+	else
+		done = 1;
+	end;    
+    it = it+1;
+%------ DEBUG
+    time_inloop = cputime-time_inloop;
+    fprintf('       time = %12.5f\n',time_inloop);
+%------ DEND
+end;
+dag	= seeddag;
+
+%----------------------------------------------------------------------
+
+outcount = 0; 
+best_score = score_dags(data,nodesizes,	{seeddag},'scoring_fn',scoring_fn);
+while outcount < 2
+	innercount = 0;
+	for	i=1:N
+		for	j=1:N
+		   if i==j,	continue;	 end;
+		   if seeddag(i,j) == 0	 % No edge i-->j, then try to add it
+			   tempdag = seeddag;
+			   tempdag(i,j)	= 1;
+			   if acyclic(tempdag)
+					temp_score = score_dags(data,nodesizes,	{tempdag},'scoring_fn',scoring_fn);
+					if temp_score >	best_score
+						seeddag	= tempdag;
+						best_score=	temp_score;
+						innercount = innercount	+1;
+					end;
+			   end
+		   else	 % exists edge i--j, then try reverse it or	remove it
+			   tempdag = seeddag;
+			   tempdag(i,j)	= 0; tempdag(j,i) =	1; 
+			   if acyclic(tempdag)
+				   temp_score =	score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+				   if temp_score > best_score
+					   seeddag = tempdag;
+					   best_score =	temp_score;
+					   innercount =	innercount +1;
+				   else
+					   tempdag = seeddag;
+					   tempdag(i,j)	= 0;
+					   temp_score =	score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+					   if temp_score > best_score
+						   seeddag = tempdag;
+						   best_score= temp_score;
+						   innercount =	innercount +1;
+					   end;
+				   end;
+			   else
+				   tempdag = seeddag;
+				   tempdag(i,j)=0;
+				   temp_score =	score_dags(data,nodesizes, {tempdag},'scoring_fn',scoring_fn);
+				   if temp_score > best_score
+					   seeddag = tempdag;
+					   best_score= temp_score;
+					   innercount =	innercount +1;
+				   end;
+			   end;
+		   end;
+		end; % end for j
+	end; % end for i
+	if innercount == 0
+		outcount = outcount	+1;
+	end;
+end;  %	end	while
+
+%------ DEBUG
+fprintf('DEBUG: Number of iterations = %d\n',it);
+fprintf('DEBUG: Outcount = %d\n',outcount);
+%------ DEND
+dag	= seeddag;
+
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabxx.m b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabxx.m
new file mode 100644
index 00000000..c0d00f05
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabxx.m
@@ -0,0 +1,141 @@
+function [dag,best_score] =	learn_struct_gs_dtabxx(data, nodesizes, seeddag, varargin)
+%
+% LEARN_STRUCT_GS(data,seeddag)	learns a structure of Bayesian net by Greedy Search.
+% dag =	learn_struct_gs(data, nodesizes, seeddag)
+%
+% dag: the final structure matrix
+% Data : training data,	data(i,m) is the m obsevation of node i
+% Nodesizes: the size array	of different nodes
+% seeddag: given seed Dag for hill climbing, optional
+%
+% by Gang Li @ Deakin University (gli73@hotmail.com)
+%
+% -----------------------------------------------------
+%
+% Modified from	learn_struct_gs (SLP 1.3)
+% to learn structure of	BN with	tabular	nodes:
+%
+% 1) make use of score decomposibility 
+% 2) replace score_family with score_family_x.c that calculates 
+%    Bayesian score with default parameters (non-adjustable)%
+% 3) replace score_dags with score_dag_x.c
+%
+% by Darima	<darrimma@yahoo.com>, 28/12/2005
+%
+
+% useold = 0;
+% if useold
+%     best_score = score_dags(data,nodesizes,	{seeddag},'scoring_fn',scoring_fn);
+% else
+%     best_score = score_dag_x(datax,nodesizes,seeddag);
+% end
+
+[N ncases] = size(data);
+if (nargin < 3 ) 
+	seeddag	= zeros(N,N); %	mk_rnd_dag(N); %call BNT function
+elseif ~acyclic(seeddag)
+	seeddag	= mk_rnd_dag(N); %zeros(N,N);
+end;
+
+% set default params (the same as in score_dags)
+for i=1:N
+  type{i} = 'tabular';
+  params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
+end
+scoring_fn = 'bayesian';
+discrete = 1:N;
+verbose	 = 'yes';
+
+% get params
+args = varargin;
+nargs =	length(args);
+if length(args)	> 0
+	if isstr(args{1})
+		for	i =	1:2:nargs
+			switch args{i}
+			case 'scoring_fn', scoring_fn =	args{i+1};
+            case 'type',       type = args{i+1}; 
+            case 'discrete',   discrete = args{i+1}; 
+            case 'params',     
+                if isempty(args{i+1}), params = cell(1,n); 
+                else params = args{i+1};  end;           
+			case 'verbose',	 verbose  =	strcmp(args{i+1},'yes');
+			end;
+		end;
+	end;
+end;
+
+%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+datax = data-ones(N,ncases);     %%+++++++++++++++++++++++ IMPORTANT!!!!
+%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+useold = 0;
+if useold
+    best_score = score_dags(data,nodesizes,	{seeddag},'scoring_fn',scoring_fn);
+else
+    best_score = score_dag_x(datax,nodesizes,seeddag);
+end
+done = 0;
+it = 0;
+
+while ~done
+	[dags,op,nodes]	= mk_nbrs_of_dag(seeddag);
+	nbrs = length(dags);
+%------ DEBUG    
+%     fprintf('DEBUG: nbrs = %6d\n',nbrs);
+%     time_inloop = cputime;
+%------ DEND
+
+%     scores	= score_dags(data, nodesizes, dags,'scoring_fn',scoring_fn);
+	scores = zeros(nbrs,1);
+	for	i =	1:nbrs
+		xj = nodes(i,2);	
+		ps_old = parents(seeddag, xj)';
+		ps_new = parents(dags{i}, xj)';
+        scor_old = score_family_x([datax(ps_old,:);datax(xj,:)],...
+            [nodesizes(ps_old),nodesizes(xj)]);
+        scor_new = score_family_x([datax(ps_new,:);datax(xj,:)],...
+            [nodesizes(ps_new),nodesizes(xj)]);
+		scores(i) =	best_score - scor_old +	scor_new;
+		if isequal(op{i},'rev')
+		    xi = nodes(i,1);
+			ps_old = parents(seeddag, xi);
+			ps_new = parents(dags{i}, xi);
+            scor_old = score_family_x([datax(ps_old,:);datax(xi,:)],...
+                [nodesizes(ps_old),nodesizes(xi)]);
+            scor_new = score_family_x([datax(ps_new,:);datax(xi,:)],...
+                [nodesizes(ps_new),nodesizes(xi)]);
+			scores(i) =	scores(i) - scor_old +	scor_new;            
+		end
+    end
+    
+	max_score =	max(scores);
+	new	= find(scores == max_score );
+
+%------ DEBUG
+%     fprintf('       -> max_score = %7.5f\n',max_score);
+%     fprintf('       -> find(scores == max_score): %d... of %d\n',...
+%         new(1),length(new));
+%------ DEND
+	if ~isempty(new) & (max_score >	best_score)
+		p =	sample_discrete(normalise(ones(1, length(new))));
+		best_score = max_score;
+		seeddag	= dags{new(p)};
+	else
+		done = 1;
+	end;
+    
+    it = it+1;
+%------ DEBUG
+%     time_inloop = cputime-time_inloop;
+%     fprintf('       time = %12.5f\n',time_inloop);
+%------ DEND
+end;
+dag	= seeddag;
+
+%------ DEBUG
+% fprintf('DEBUG: Number of iterations = %d\n',it);
+% fprintf('DEBUG: Outcount = %d\n',0);
+%------ DEND
+
+
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_dag_x.c b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_dag_x.c
new file mode 100644
index 00000000..d72bf591
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_dag_x.c
@@ -0,0 +1,72 @@
+//
+// File: score_dag_x.c
+//
+// MATLAB:
+// function score = score_dag_x(data,sz,dag)
+//
+// DESCRIPTION:
+// Calculates Bayesian score for the DAG with tabular nodes
+// (as in log_marg_prob_node.m):
+//      data        [nsz x ndata] array
+//                  IMPORTANT:                    <--- !!!!
+//                  - parents go first                             
+//                  - values of the nodes are in the range 0..sz[i] 
+//      sz          sizes of the nodes
+//		dag			DAG
+//      ndata       # of observations
+//      nsz         # of nodes in the DAG
+// with default parameters:  
+//      params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }
+//                   ...and 'dirichlet_type','BDeu'
+//      scoring_fn = 'bayesian';
+//      etc...
+//
+// EXAMPLE:
+//   score_family_x([1 2 1; 1 2 2; 1 1 1],[2 2 2]);
+//
+// by Darima <darrimma@yahoo.com>, 27/12/2005 
+//
+
+#include "mex.h"
+#include "score_x.h"
+
+#define	IN_DATA	    prhs[0]
+#define IN_SZ       prhs[1]
+#define IN_DAG      prhs[2]
+#define	OUT_SCORE	plhs[0]
+#if !defined(MAX)
+#define	MAX(A, B)	((A) > (B) ? (A) : (B))
+#endif
+
+void mexFunction( int nlhs, mxArray *plhs[],
+                  int nrhs, const mxArray *prhs[] )   
+{ 
+    int nsz, ndata;
+    double *sz, *data,*dag, *score;
+   
+    /* Check for proper number of arguments */
+    
+    if (nrhs != 3) { 
+        mexErrMsgTxt("Three input arguments required."); 
+    } else if (nlhs > 1) {
+        mexErrMsgTxt("Too many output arguments."); 
+    } 
+
+    /* Assign input arguments */ 
+    
+    data   = mxGetPr(IN_DATA);
+    sz     = mxGetPr(IN_SZ);
+	dag	   = mxGetPr(IN_DAG);
+	ndata  = mxGetN(IN_DATA); 
+    nsz    = MAX(mxGetM(IN_SZ),mxGetN(IN_SZ)); 
+	
+    /* Create return argument */ 
+    
+    OUT_SCORE = mxCreateDoubleScalar(mxREAL);
+    score = mxGetPr(OUT_SCORE);
+
+    /* Do the actual computations in a subroutine */
+    
+    score_dag_x(data,sz,dag,ndata,nsz,score);
+    return;
+}
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_family_x.c b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_family_x.c
new file mode 100644
index 00000000..ce45fd94
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_family_x.c
@@ -0,0 +1,71 @@
+//
+// File: score_family_x.c
+//
+// MATLAB:
+// function score = score_family_x(data,sz)
+//
+// DESCRIPTION:
+// Calculates Bayesian score for the family of tabular nodes
+// (as in log_marg_prob_node.m):
+//      data        [nsz x ndata] array
+//                  IMPORTANT:                    <--- !!!!
+//                  - parents go first                             
+//                  - values of the nodes are in the range 0..sz[i] 
+//      sz          sizes of the nodes in the family
+//      ndata       # of observations
+//      nsz         # of nodes in the family
+// with default parameters:  
+//      params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }
+//                   ...and 'dirichlet_type','BDeu'
+//      scoring_fn = 'bayesian';
+//      etc...
+//
+// EXAMPLE:
+//   score_family_x([1 2 1; 1 2 2; 1 1 1],[2 2 2]);  // ans = -2.0794
+//
+// by Darima <darrimma@yahoo.com>, 27/12/2005 
+//
+
+#include "mex.h"
+#include "score_x.h"
+
+#define	IN_DATA	    prhs[0]
+#define IN_SZ       prhs[1]
+#define	OUT_SCORE	plhs[0]
+#if !defined(MAX)
+#define	MAX(A, B)	((A) > (B) ? (A) : (B))
+#endif
+
+void mexFunction( int nlhs, mxArray *plhs[],
+                  int nrhs, const mxArray *prhs[] )   
+{ 
+    int nsz, ndata;
+    double  *sz, *data, *score;
+	int i,j;
+    
+    /* Check for proper number of arguments */
+    
+    if (nrhs != 2) { 
+        mexErrMsgTxt("Two input arguments required."); 
+    } else if (nlhs > 1) {
+        mexErrMsgTxt("Too many output arguments."); 
+    } 
+
+    /* Assign input arguments */ 
+    
+    data   = mxGetPr(IN_DATA);
+    sz     = mxGetPr(IN_SZ);
+	ndata  = mxGetN(IN_DATA); 
+    nsz    = MAX(mxGetM(IN_SZ),mxGetN(IN_SZ));   
+
+    /* Create return argument */ 
+    
+    OUT_SCORE = mxCreateDoubleScalar(mxREAL);
+    score = mxGetPr(OUT_SCORE);
+
+    /* Do the actual computations in a subroutine */
+    
+    score_family_x(data,sz,nsz,ndata,score);
+    return;
+}
+
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.c b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.c
new file mode 100644
index 00000000..b63111f6
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.c
@@ -0,0 +1,168 @@
+//
+// File: score_x.c
+//
+//
+// by Darima <darrimma@yahoo.com>, 27/12/2005 
+//
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <malloc.h>
+#include "score_x.h"
+
+void score_dag_x(double *data,double *sz,double *dag,
+    int ndata,int nsz,double* score)
+//
+// GOTO: score_dag_x.c
+// by Darima <darrimma@yahoo.com>, 27/12/2005 
+//    
+{
+	double *fdata,*fsz,*fscore;
+    int fnsz,*family; 
+	int i,j,k; 
+        
+	fdata = (double*)malloc(nsz*ndata*sizeof(double));
+	fsz = (double*)malloc(nsz*sizeof(double));
+    fscore = (double*)malloc(sizeof(double));
+    family = (int*)malloc(nsz*sizeof(int));
+	*score = 0;      
+   
+	for (j=0;j<nsz;j++) {
+        
+		// find parents of current node
+		        
+        fnsz = 0;
+		for (i=0;i<nsz;i++)
+			if (dag[i+j*nsz]==1) {family[fnsz] = i; fnsz++;}
+        family[fnsz] = j; fnsz++;
+        
+        // initialize data for current family
+        
+		for (i=0;i<fnsz;i++)
+            for (k=0;k<ndata;k++) {
+                fdata[i+k*fnsz] = data[family[i]+k*nsz];
+                fsz[i] = sz[family[i]];
+            }
+        
+		// calculater score for the family
+        
+        score_family_x(fdata,fsz,fnsz,ndata,fscore);
+        *score += *fscore;
+	}
+
+	free(fdata); free(fsz); free(fscore); free(family);
+	return;
+}
+
+void score_family_x(double *data,double *sz,
+    int nsz,int ndata,double* score)
+//
+// GOTO: score_family_x.c
+// by Darima <darrimma@yahoo.com>, 27/12/2005 
+//
+{ 
+    int self_sz, ps_sz;
+    double prior1, prior2, *count;
+    int i,j,k; double idx,tsz,N_ij;
+        
+    // self_sz - number of values of current node
+    // ps_sz   - number of configurations of parents    
+    
+    self_sz = sz[nsz-1]; 
+    ps_sz = 1; for (i=0;i<nsz-1;i++){ps_sz *= sz[i];}
+       
+    // calculate counts
+    //
+    // example: two binary parents for binary node
+    //  idx |  p1 p2 |   0   1   <-- values of the node
+    //  ---------------------
+    //   0  |  0  0  |   1   0
+    //   1  |  1  0  |   0  83
+    //   2  |  0  1  |   3   0   <-- counts
+    //   3  |  1  1  |  13   0
+    //  ^------------------------ configurations of parents
+    // count = [ps_sz x self_sz] array
+        
+    count = (double*)malloc(ps_sz*self_sz*sizeof(double));
+    for (i=0;i<self_sz*ps_sz;i++) { count[i] = 0; }
+
+    for (j=0;j<ndata;j++) {
+        
+        // for every data case:
+        // 1) calculate parent configuration index: p1_val,p2_val -> idx
+        // 2) increment corresponding count: count[idx,node_val]++
+        //
+        // example:
+        //             case1  case2  etc
+        // -----------------------------
+        // p1   |      0      1      ...
+        // p2   |      0      0      ...
+        // node |      1      1      ...
+        // 
+        // notice: array[i,j,k] -> array[i+j*ni+k*ni*nj]
+        //
+        // example 1: data[i,j] -> data[i+j*nsz]
+        //
+        // example 2: consider family with 5 parents
+        // - [p1,p2,..,p5] -> IDX5 
+        //   IDX5 = p1 + p2*p1_sz + p3*p1_sz*p2_sz + ... +
+        //          p5*p1_sz*p2_sz*p3_sz*p4_sz
+        // - [p1,p2,..,p5,node_val] -> IDXFAM
+        //   IDXFAM = IDX5 + node_val*ps_sz
+        //   count[IDX5,node_val] ->count[IDXFAM]
+        
+        idx = data[0+j*nsz]; tsz = 1;
+        for (i=1;i<nsz;i++) {tsz *= sz[i-1]; idx += data[i+j*nsz]*tsz;}
+        count[(int)idx]++;
+    }
+       
+    // calculate priors (BDeu,1)
+    
+    prior1 = 1/(double)(self_sz*ps_sz);
+    prior2 = 1/(double)(ps_sz);
+       
+    // calculate score
+    //
+    // LL = log[  PROD_j gamma(alpha_ij)/gamma(alpha_ij + N_ij)  *
+	//            PROD_k gamma(alpha_ijk + N_ijk)/gamma(alpha_ijk)  ] =
+    //    = SUM_j {  gammaln(alpha_ij)-gammaln(alpha_ij + N_ij) +
+    //         SUM_k [ gammaln(alpha_ijk + N_ijk)-gammaln(alpha_ijk) ]  }
+    //
+    // alpha_ijk  <- prior1
+    // alpha_ij   <- prior2
+    // N_ijk      <- count[j,k]
+
+    *score = 0;
+    for (j=0;j<ps_sz;j++) {    
+        N_ij = .0; 
+        for (k=0;k<self_sz;k++) {
+            N_ij += count[j+k*ps_sz];
+            *score += gammaln(prior1+count[j+k*ps_sz])-gammaln(prior1);
+        }
+        *score += gammaln(prior2)-gammaln(prior2 + N_ij);
+    }
+
+    free(count);
+    return;
+}
+
+double gammaln(double xx)
+//
+// Returns the value ln[Gamma(xx)] for xx > 0.
+// from "Numerical Recipes in C"
+//
+{
+    double x,y,tmp,ser;
+    static double cof[6]={
+        76.18009172947146,-86.50532032941677,
+        24.01409824083091,-1.231739572450155,
+        0.1208650973866179e-2,-0.5395239384953e-5};
+    int j;
+    y=x=xx;
+    tmp=x+5.5;
+    tmp -= (x+0.5)*log(tmp);
+    ser=1.000000000190015;
+    for (j=0;j<=5;j++) ser += cof[j]/++y;
+    return -tmp+log(2.5066282746310005*ser/x);
+}
diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.h b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.h
new file mode 100644
index 00000000..8fc1c213
--- /dev/null
+++ b/sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.h
@@ -0,0 +1,19 @@
+//
+// File: score_x.h
+//
+//
+// by Darima <darrimma@yahoo.com>, 27/12/2005 
+//
+
+#ifndef score_x_h
+#define mex_h
+
+void score_family_x(double *data,double *sz,
+    int ndata,int nsz,double* score);
+
+void score_dag_x(double *data,double *sz,double *dag,
+    int nsz,int ndata,double* score);
+
+double gammaln(double xx);
+
+#endif // score_x_h