From 8070dc963753142bb86c4ed698d91fd623ed28e7 Mon Sep 17 00:00:00 2001 From: ziejd2 Date: Thu, 28 Sep 2017 15:04:40 -0500 Subject: 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 --- .../scoring/calculate_mutual_information_array.m | 32 +++ .../bnt-master/SLP/scoring/cond_indep_chisquare.m | 184 +++++++++++++++++ .../SLP/scoring/cond_mutual_info_score.m | 28 +++ sourcecodes/bnt-master/SLP/scoring/kl_divergence.m | 48 +++++ .../bnt-master/SLP/scoring/kl_divergence2.m | 43 ++++ .../bnt-master/SLP/scoring/mutual_info_score.m | 25 +++ .../bnt-master/SLP/scoring/score_add_to_cache.m | 67 +++++++ .../SLP/scoring/score_dag.c/learn_struct_gs_dtab.m | 179 +++++++++++++++++ .../score_dag.c/learn_struct_gs_dtab_INFO.txt | 91 +++++++++ .../scoring/score_dag.c/learn_struct_gs_dtabx.m | 181 +++++++++++++++++ .../scoring/score_dag.c/learn_struct_gs_dtabxx.m | 141 +++++++++++++ .../SLP/scoring/score_dag.c/score_dag_x.c | 72 +++++++ .../SLP/scoring/score_dag.c/score_family_x.c | 71 +++++++ .../bnt-master/SLP/scoring/score_dag.c/score_x.c | 168 ++++++++++++++++ .../bnt-master/SLP/scoring/score_dag.c/score_x.h | 19 ++ sourcecodes/bnt-master/SLP/scoring/score_dags.m | 74 +++++++ sourcecodes/bnt-master/SLP/scoring/score_family.m | 218 +++++++++++++++++++++ .../bnt-master/SLP/scoring/score_find_in_cache.m | 63 ++++++ .../bnt-master/SLP/scoring/score_init_cache.m | 29 +++ 19 files changed, 1733 insertions(+) create mode 100644 sourcecodes/bnt-master/SLP/scoring/calculate_mutual_information_array.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/cond_indep_chisquare.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/cond_mutual_info_score.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/kl_divergence.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/kl_divergence2.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/mutual_info_score.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_add_to_cache.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtab_INFO.txt create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabx.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/learn_struct_gs_dtabxx.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_dag_x.c create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_family_x.c create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.c create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dag.c/score_x.h create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_dags.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_family.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_find_in_cache.m create mode 100644 sourcecodes/bnt-master/SLP/scoring/score_init_cache.m (limited to 'sourcecodes/bnt-master/SLP/scoring') diff --git a/sourcecodes/bnt-master/SLP/scoring/calculate_mutual_information_array.m b/sourcecodes/bnt-master/SLP/scoring/calculate_mutual_information_array.m new file mode 100644 index 00000000..1128e411 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/calculate_mutual_information_array.m @@ -0,0 +1,32 @@ +function [mi] = calculate_mutual_information_array(data) +% FUNCTION [MI_ARRAY] = CALCULATE_MUTUAL_INFORMATION_ARRAY(DATA) +% calculates the mutual information between all pairs of variables +% Data must be discrete, and take values 1,2,...,size +% data(i,m) is the node i in the case m. + +[num_nodes num_examples] = size(data); + +node_sizes = max(data'); +for i = 1:num_nodes + for ic = 1:node_sizes(i) % I CLASS ic + px(i,ic) = sum(data(i,:)==ic); + for j = 1:num_nodes % J CLASS jc + for jc = 1:node_sizes(j) + pxy(i,ic,j,jc) = sum( (data(i,:)==ic) & (data(j,:)==jc) ); + end + mi(i,j) = 0; + end + end +end + +for i = 1:num_nodes + for ic = 1:node_sizes(i) + for j = 1:num_nodes + for jc = 1:node_sizes(j) + if( pxy(i,ic,j,jc)~=0 & px(i,ic)~=0 & px(j,jc)~= 0) + mi(i,j) = mi(i,j) + pxy(i,ic,j,jc)*log2( num_examples*pxy(i,ic,j,jc)/(px(i,ic)*px(j,jc)) )/num_examples; + end + end + end + end +end diff --git a/sourcecodes/bnt-master/SLP/scoring/cond_indep_chisquare.m b/sourcecodes/bnt-master/SLP/scoring/cond_indep_chisquare.m new file mode 100644 index 00000000..efc162cc --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/cond_indep_chisquare.m @@ -0,0 +1,184 @@ +function [CI, Chi2, alpha2] = cond_indep_chisquare(X, Y, S, Data, test, alpha, ns) +% COND_INDEP_CHISQUARE Test if X indep Y given Z +% using either chisquare test or likelihood ratio test G2 +% +% [CI Chi2 Prob_Chi2] = cond_indep_chisquare(X, Y, S, Data, test, alpha, node_sizes) +% +% Input : +% Data is the data matrix, N columns * NbVar rows +% X is the index of variable X in Data matrix +% Y is the index of variable Y in Data matrix +% S are the indexes of variables in set S +% alpha is the significance level (default: 0.05) +% test = 'pearson' for Pearson's chi2 test +% 'LRT' for G2 likelihood ratio test (default) +% node_sizes (default: max(Data')) +% +% Output : +% CI = test result (1=conditional independency, 0=no) +% Chi2 = chi2 value (-1 if not enough data to perform the test --> CI=0) +% +% +% V1.4 : 24 july 2003 (Ph. Leray - philippe.leray@univ-nantes.fr) +% +% +% Things to do : +% - do not use 'find' in nij computation (when S=empty set) +% - find a better way than 'warning off/on' in tmpij, tmpijk computation +% + +if nargin < 5, test = 'LRT'; end +if nargin < 6, alpha = 0.05; end +if nargin < 7, ns = max(Data'); end + +Data=Data'; + +N = size(Data,1); +qi=ns(S); +tmp=[1 cumprod(qi(1:end-1))]; +qs=1+(qi-1)*tmp'; +if isempty(qs), + nij=zeros(ns(X),ns(Y)); + df=prod(ns([X Y])-1)*prod(ns(S)); +else + +% Commented by Mingyi +% nijk=zeros(ns(X),ns(Y),qs); +% tijk=zeros(ns(X),ns(Y),qs); +% Commention ends +% Added by Mingyi + nijk=zeros(ns(X),ns(Y),1); + tijk=zeros(ns(X),ns(Y),1); +% Addition ends + df=prod(ns([X Y])-1)*qs; +end + + +if (N<10*df) + % Not enough data to perform the test + Chi2=-1; + CI=0; + +elseif isempty(S) + for i=1:ns(X), + for j=1:ns(Y), + nij(i,j)=length(find((Data(:,X)==i)&(Data(:,Y)==j))) ; + end + end + restr=find(sum(nij,1)==0); + if ~isempty(restr) + nij=nij(:,find(sum(nij,1))); + end + + tij=sum(nij,2)*sum(nij,1)/N ; + + switch test + case 'pearson', + tmpij=nij-tij; + + [xi yj]=find(tij<10); + for i=1:length(xi), + tmpij(xi(i),yj(i))=abs(tmpij(xi(i),yj(i)))-0.5; + end + + warning off; + tmp=(tmpij.^2)./tij; + warning on; + tmp(find(tmp==Inf))=0; + + case 'LRT', + warning off; + tmp=nij./tij; + warning on; + tmp(find(tmp==Inf | tmp==0))=1; + tmp(find(tmp~=tmp))=1; + tmp=2*nij.*log(tmp); + + otherwise, + error(['unrecognized test ' test]); + end + + Chi2=sum(sum(tmp)); + alpha2=1-chisquared_prob(Chi2,df); + CI=(alpha2>=alpha) ; + +else + SizeofSSi=1; + for exemple=1:N, + i=Data(exemple,X); + j=Data(exemple,Y); + Si=Data(exemple,S)-1; + %Added by Mingyi + if exemple==1 + SSi(SizeofSSi,:)=Si; + nijk(i,j,SizeofSSi)=1; + else + flag=0; + for iii=1:SizeofSSi + if isequal(SSi(iii,:),Si) + nijk(i,j,iii)=nijk(i,j,iii)+1; + flag=1; + end + end + if flag==0 + SizeofSSi=SizeofSSi+1; + SSi(SizeofSSi,:)=Si; + nijk(i,j,SizeofSSi)=1; + end + end + %Addition ends + %Commented by Mingyi +% k=1+Si*tmp'; +% nijk(i,j,k)=nijk(i,j,k)+1; + %Commention ends + end + + nik=sum(nijk,2); + njk=sum(nijk,1); + N2=sum(njk); + + % for k=1:qs, %Commented by Mingyi + for k=1:SizeofSSi %Added by Mingyi + if N2(:,:,k)==0 + tijk(:,:,k)=0; + else + tijk(:,:,k)=nik(:,:,k)*njk(:,:,k)/N2(:,:,k); + end + end + + switch test + case 'pearson', + tmpijk=nijk-tijk; + + [xi yj]=find(tijk<10); + for i=1:length(xi), + tmpijk(xi(i),yj(i))=abs(tmpijk(xi(i),yj(i)))-0.5; + end + + warning off; + tmp=(tmpijk.^2)./tijk; + warning on; + tmp(find(tmp==Inf))=0; + + case 'LRT', + warning off; + tmp=nijk./tijk; + warning on; + tmp(find(tmp==Inf | tmp==0))=1; + tmp(find(tmp~=tmp))=1; + tmp=2*nijk.*log(tmp); + + otherwise, + error(['unrecognized test ' test]); + end + + Chi2=sum(sum(sum(tmp))); + alpha2=1-chisquared_prob(Chi2,df); + CI=(alpha2>=alpha) ; + +end +clear tijk +clear nijk +clear nij +clear tij +clear tmpijk \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/scoring/cond_mutual_info_score.m b/sourcecodes/bnt-master/SLP/scoring/cond_mutual_info_score.m new file mode 100644 index 00000000..b43cbfeb --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/cond_mutual_info_score.m @@ -0,0 +1,28 @@ +function score = cond_mutual_info_score(i,si,j,sj,c,sc,data) +% G = cond_mutual_info_score(i,si,j,sj,c,sc,data) +% Only for tabular node which values are 1,2,...,size . +% si is size of node i, sj is size of node j, sc is the size of node c. +% data(i,m) is the node i in the case m. +% +% +% pphilippe.leray@univ-nantes.fr, francois.olivier.c.h@gmail.com + +[n N]=size(data); +Pc=hist(data(c,:),1:sc)/N; +score=0; + +for cvalue=1:sc, + ind=find(data(c,:)==cvalue); + Nj=hist(data(j,ind),1:sj); + Ni=hist(data(i,ind),1:si); + NiNj=Ni'*Nj; + + for k=1:si + ind2=find(data(i,ind)==k) ; + Nij(k,:) = hist(data(j,ind(ind2)),1:sj); + end + + % sommons les valeurs non-infinies: + ind=find(NiNj~=0 & Nij~=0); + score=score+Pc(cvalue)*sum(sum(Nij(ind).*log(N*Nij(ind)./NiNj(ind))/N)); +end diff --git a/sourcecodes/bnt-master/SLP/scoring/kl_divergence.m b/sourcecodes/bnt-master/SLP/scoring/kl_divergence.m new file mode 100644 index 00000000..f4e18147 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/kl_divergence.m @@ -0,0 +1,48 @@ +function KLdiv = KL_divergence(bnetP, bnetQ) +% KL_DIVERGENCE computes the Kullback-Leibler divergence between two BNET distributions +% KLdiv = KL_divergence(bnetP, bnetQ) +% +% Output : +% div = sum_x P(x).log(P(x)/Q(x)) +% +% Rem : +% This version is optimized for speed, but can use too many memory +% ==> if you have a memory problem, use kl_divergence2 instead +% + +% ONLY FOR TABULAR NODES +% Make sure that you have done the params learning. +% +% V1.1 : 8 oct 2004 (Ph. Leray - philippe.leray@univ-nantes.fr) + +N = size(bnetP.dag,1); +N2 = size(bnetQ.dag,1); +ns= bnetP.node_sizes; +ns2= bnetQ.node_sizes; +if N~=N2, error('size of dags must be the same'), end +if ns~=ns2, error('node sizes of dags must be the same'), end +tiny = exp(-700); +KLdiv=0; + +inst = ind2subv(ns, 1:prod(ns)); + %Px=1; Qx=1; + for i=1:N, + ps = parents(bnetP.dag, i); + %e = bnetP.equiv_class(i); + %[tmp Px(:,i)] = prob_node(bnetP.CPD{e}, inst(:,i)', inst(:,ps)'); + [tmp Px(:,i)] = prob_node(bnetP.CPD{i}, inst(:,i)', inst(:,ps)'); + + ps = parents(bnetQ.dag, i); + %e = bnetQ.equiv_class(i); + %[tmp Qx(:,i)] = prob_node(bnetQ.CPD{e}, inst(:,i)', inst(:,ps)'); + [tmp Qx(:,i)] = prob_node(bnetQ.CPD{i}, inst(:,i)', inst(:,ps)'); + end + Px=prod(Px,2); + Px = Px + (Px==0)*tiny; % replace 0s by tiny + Qx=prod(Qx,2); + Qx = Qx + (Qx==0)*tiny; % replace 0s by tiny + + %%%%% Faut-il diviser par le nb de configurations possibles ? (sum => mean) + KLdiv = sum(Px.*log(Px./Qx)); + %end + diff --git a/sourcecodes/bnt-master/SLP/scoring/kl_divergence2.m b/sourcecodes/bnt-master/SLP/scoring/kl_divergence2.m new file mode 100644 index 00000000..7b9664ba --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/kl_divergence2.m @@ -0,0 +1,43 @@ +function KLdiv = KL_divergence2(bnetP, bnetQ) +% KL_DIVERGENCE2 computes the Kullback-Leibler divergence between two BNET distributions +% KLdiv = KL_divergence2(bnetP, bnetQ) +% +% Output : +% div = sum_x P(x).log(P(x)/Q(x)) +% +% Rem : +% This version is optimized for memory use, but quite slow !!! +% ==> if you have no memory problem, use kl_divergence instead +% +% ONLY FOR TABULAR NODES +% Make sure that you have done the params learning. +% +% V1.1 : 8 oct 2004 (Ph. Leray - philippe.leray@univ-nantes.fr) + +N = size(bnetP.dag,1); +N2 = size(bnetQ.dag,1); +ns= bnetP.node_sizes; +ns2= bnetQ.node_sizes; +if N~=N2, error('size of dags must be the same'), end +if ns~=ns2, error('node sizes of dags must be the same'), end +tiny = exp(-700); +KLdiv=0; + +for i=1:prod(ns), + inst = ind2subv(ns, i); % i'th instantiation + Px=1; Qx=1; + for i=1:N, + ps = parents(bnetP.dag, i); + e = bnetP.equiv_class(i); + [tmp Pxi] = prob_node(bnetP.CPD{e}, inst(i), inst(ps)'); + Px=Px*Pxi; + ps = parents(bnetQ.dag, i); + e = bnetQ.equiv_class(i); + [tmp Qxi] = prob_node(bnetQ.CPD{e}, inst(i), inst(ps)'); + Qx=Qx*Qxi; + end + Px = Px + (Px==0)*tiny; % replace 0s by tiny + Qx = Qx + (Qx==0)*tiny; % replace 0s by tiny + KLdiv = KLdiv + Px*log(Px/Qx); +end + diff --git a/sourcecodes/bnt-master/SLP/scoring/mutual_info_score.m b/sourcecodes/bnt-master/SLP/scoring/mutual_info_score.m new file mode 100644 index 00000000..8808a124 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/mutual_info_score.m @@ -0,0 +1,25 @@ +function score = mutual_info_score(i,si,j,sj,data) +% G = mutual_info_score(i,si,j,sj,data) +% Only for tabular node which values are 1,2,...,size . +% si is size of node i, sj is size of node j. +% data(i,m) is the node i in the case m. +% +% Ref : +% C. Chow and C. Liu (1968). Approximating discrete probability distributions with dependence trees. +% IEEE Transactions on Information Theory, 14(3):462--467, May 1968. +% +% francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr, wangxiangyang@sjtu.edu.cn + +[n N]=size(data); +Nj=hist(data(j,:),1:sj); +Ni=hist(data(i,:),1:si); +NiNj=Ni'*Nj; + +for k=1:si + ind=find(data(i,:)==k) ; + Nij(k,:) = hist(data(j,ind),1:sj); +end + +% sommons les valeurs non-infinies: +ind=find(NiNj~=0 & Nij~=0); +score=sum(sum(Nij(ind).*log(N*Nij(ind)./NiNj(ind))/N)); diff --git a/sourcecodes/bnt-master/SLP/scoring/score_add_to_cache.m b/sourcecodes/bnt-master/SLP/scoring/score_add_to_cache.m new file mode 100644 index 00000000..3c68e4d8 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/score_add_to_cache.m @@ -0,0 +1,67 @@ +function [cache, place] = score_add_to_cache(cache,j,ps,score,scoring_fn) +% [cache place] = score_add_to_cache(cache,j,ps,score,scoring_fn) +% +% j is the son node, +% ps is the list of parents of j, for example [12 5 7], +% score is the score to add for this familly. +% scoring_fn is 'bic' or 'bayesian'. +% +% place = where the entry was add. +% +% example for 2 nodes with cache of size 5 : +% +% cache = +% 5 b 0 0 0 --> number of writing in cache +1 and b==1 iff the cache is full +% 0 0 1 -239.12 1 --> 1st familly in the cache (node 1 without parents) calculate with bic +% 0 0 2 -318.98 1 +% 1 0 2 -189.23 2 --> 3rd familly in the cache (node 2 with 1 as parent) calculate with bayesian +% 0 1 1 -251.09 1 +% .ps2bool. j score 1or2 --> new entry +% | | | | | +% | | | | |___> 1 for 'bic' or 2 for 'bayesian' +% | | | |__________> score of the familly +% | | |_________________> son node of the familly +% | |__________________________> ==1 iff node 2 is parent of son node +% |______________________________> ==1 iff node 1 is parent of son node +% +% If the cache is FULL then the new place is RanDoMly choose. +% +% francois.olivier.c.h@gmail.com + +N=size(cache,2)-3; +place=0; + +if ~isempty(find(ps==j)) + disp('This is a cyclic entry, nothing was done.'); +elseif j>N | j<0 + disp('This entry is not valid, nothing was done.'); +else + + switch scoring_fn + case 'bic', + fn=1; + case 'bayesian', + fn=2; + otherwise, + fn=3; + %error(['unrecognized scoring fn ' scoring_fn]); + end + L=size(cache,1); + + if cache(1,2)==0 + place=cache(1,1); + else + place=ceil(rand(1)*(L-1))+1; + end + + cache(place,:)=0 + cache(place,ps)=1; + cache(place,N+1)=j; + cache(place,N+2)=score; + cache(place,N+3)=fn; + cache(1,1)=place+1; + if place==L | cache(1,2)~=0 + cache(1,2)=1 + end + +end 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 , 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" +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 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 , 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 , 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 , 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 , 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 , 27/12/2005 +// + +#include +#include +#include +#include +#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 , 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, 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 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 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 , 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 diff --git a/sourcecodes/bnt-master/SLP/scoring/score_dags.m b/sourcecodes/bnt-master/SLP/scoring/score_dags.m new file mode 100644 index 00000000..acda3e0f --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/score_dags.m @@ -0,0 +1,74 @@ +function [score, cache] = score_dags(data, ns, dags, varargin) +% SCORE_DAGS Compute the score of one or more DAGs +% score = score_dags(data, ns, dags, varargin) +% +% data{i,m} = value of node i in case m (can be a cell array). +% node_sizes(i) is the number of size of node i. +% dags{g} is the g'th dag +% score(g) is the score of the i'th dag +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ] +% Currently, only networks with all tabular nodes support Bayesian scoring. +% type - type{i} is the type of CPD to use for node i, where the type is a string +% of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ] +% params - params{i} contains optional arguments passed to the CPD constructor for node i, +% or [] if none. [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ] +% discrete - the list of discrete nodes [ 1:N ] +% clamped - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ] +% cache - data structure used to memorize local score computations (cf. SCORE_INIT_CACHE function) [ [] ] +% +% e.g., score = score_dags(data, ns, mk_all_dags(n), 'scoring_fn', 'bic', 'params', [],'cache',cache); +% +% (Caching implementation : francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr) +% ("Clamped" optimisation : Derek Hoiem ) + +[n ncases] = size(data); + +% set default params +type = cell(1,n); +params = cell(1,n); +cache=[]; +for i=1:n + type{i} = 'tabular'; + params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 }; +end +scoring_fn = 'bayesian'; +discrete = 1:n; + +isclamped = 0; % DWH +clamped = zeros(n, ncases); +u = [1:ncases]'; % DWH + +args = varargin; +nargs = length(args); +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 'clamped', clamped = args{i+1}, isclamped = 1; %DWH + case 'params', if isempty(args{i+1}), params = cell(1,n); else params = args{i+1}; end + case 'cache', cache=args{i+1} ; + end +end + +NG = length(dags); +score = zeros(1, NG); + +for j=1:n + if isclamped %DWH + u = find(clamped(j,:)==0); + end + for g=1:NG + if isempty(dags{g}) + score(g)=-Inf; + else + ps = parents(dags{g}, j); + [scor cache] = score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j}, cache); + score(g) = score(g) + scor; + end + end +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/scoring/score_family.m b/sourcecodes/bnt-master/SLP/scoring/score_family.m new file mode 100644 index 00000000..085fcf48 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/score_family.m @@ -0,0 +1,218 @@ +function [score, cache] = score_family(j, ps, node_type, scoring_fn, ns, discrete, data, args, cache) +% SCORE_FAMILY_COMPLETE Compute the score of a node and its parents given completely observed data +% score = score_family(j, ps, node_type, scoring_fn, ns, discrete, data, args, cache) +% +% data(i,m) is the value of node i in case m (can be a cell array, if contain missing value, it uses only available complete cases for an entry) +% args is a cell array containing optional arguments passed to the constructor, +% or is [] if none +% cache is a data structure used to memorize local score computations +% (cf. SCORE_INIT_CACHE function) +% +% We create a whole Bayes net which only connects parents to node, +% where node has a CPD of the specified type (with default parameters). +% We then evaluate its score ('bic' or 'bayesian') +% We should use a cache to avoid unnecessary computation. +% In particular, log_marginal_prob_node for tabular CPDs calls gammaln +% and compute_counts, both of which are slow. +% +% (Caching implementation : ofrancois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr) +% + +if (nargin<9 || isempty(cache)) , c=0; else c=1; end +%tic +if c==1 + [b,score]=score_find_in_cache(cache,j,ps,scoring_fn); +else + b=0; +end +%Tfind=toc +ccc = iscell(data); +ps = unique(ps); + +if b==0 + misv = -9999; + if ccc, data = bnt_to_mat(data,misv); end + [n ncases] = size(data); + dag = zeros(n,n); + % SML added to sort ps b/c mk_bnet, learn_params use sorted ps to make + % CPTs % Kevin had: if ~isempty(ps), dag(ps, j) = 1; end + if ~isempty(ps), dag(ps, j) = 1;, ps = sort(ps);, end + + bnet = mk_bnet(dag, ns, 'discrete', discrete); + fname = sprintf('%s_CPD', node_type); + if isempty(args) + bnet.CPD{j} = feval(fname, bnet, j); + else + bnet.CPD{j} = feval(fname, bnet, j, args{:}); + end + %tic + switch scoring_fn + case 'bic', + fam = [ps j]; + if ccc, + [tmp, available_case] = find(data(fam,:)==misv); + available_case = mysetdiff(1:ncases, available_case); + else available_case = 1:ncases; + end + bnet.CPD{j} = learn_params(bnet.CPD{j}, fam, data(:,available_case), ns, bnet.cnodes); + %L = log_prob_node(bnet.CPD{j}, data(j,:), data(ps,:)); + L = log_prob_node(bnet.CPD{j}, data(j,available_case), data(ps,available_case)); + S = struct(bnet.CPD{j}); % violate object privacy + score = L - 0.5*S.nparams*log(length(available_case)); + case 'bicmod', + fam = [ps j]; + if ccc, + [tmp, available_case] = find(data(fam,:)==misv); + available_case = mysetdiff(1:ncases, available_case); + else available_case = 1:ncases; + end + bnet.CPD{j} = learn_params(bnet.CPD{j}, fam, data(:,available_case), ns, bnet.cnodes); + L = log_prob_node(bnet.CPD{j}, data(j,available_case), data(ps,available_case)); + S = struct(bnet.CPD{j}); % violate object privacy + score = L - S.nparams*log(length(available_case)); + case 'bayesian', + fam = [ps j]; + if ccc, + [tmp, available_case] = find(data(fam,:)==misv); + available_case = mysetdiff(1:ncases, available_case); + else available_case = 1:ncases; + end + score = log_marg_prob_node(bnet.CPD{j}, data(j,available_case), data(ps,available_case)); + otherwise, + error(['unrecognized scoring fn ' scoring_fn]); + end + %Tcalc=toc + %tic + if c==1 +% fprintf('a\n') + cache=score_add_to_cache(cache,j,ps,score,scoring_fn); + end + %Tecr=toc +% else +% fprintf('*\n') +end + +%===========================Inner functions + +function [cache, place] = score_add_to_cache(cache,j,ps,score,scoring_fn) +% [cache place] = score_add_to_cache(cache,j,ps,score,scoring_fn) +% +% j is the son node, +% ps is the list of parents of j, for example [12 5 7], +% score is the score to add for this familly. +% scoring_fn is 'bic' or 'bayesian'. +% +% place = where the entry was add. +% +% example for 2 nodes with cache of size 5 : +% +% cache = +% Nw b 0 0 0 --> Nw=number of writings in cache (+1) and b==1 iff the cache is full +% 0 0 1 -239.12 1 --> 1st familly in the cache (node 1 without parents) calculate with bic +% 0 0 2 -318.98 1 +% 1 0 2 -189.23 2 --> 3rd familly in the cache (node 2 with 1 as parent) calculate with bayésian +% 0 1 1 -251.09 1 +% .ps2bool. j score 1or2 --> new entry +% | | | | | +% | | | | |___> scoring function : 1 for 'bic' or 2 for 'bayesian' +% | | | |__________> score of the familly +% | | |_________________> son node of the familly +% | |__________________________> ==1 iff node 2 is parent of son node +% |______________________________> ==1 iff node 1 is parent of son node +% +% If the cache is FULL then the new place is RanDomly choose. +% +% V1.1 : 5 may 2003 (O. Francois, Ph. Leray) + +N=size(cache,2)-3; +L=size(cache,1)-1; +cache_full=cache(1,2) ; + +place=0; + +if ismember(j,ps) + disp('This is a cyclic entry, nothing was done.'); +elseif j>N || j<=0 + disp('This entry is not valid, nothing was done.'); +else + switch scoring_fn + case 'bic', + fn=1; + case 'bayesian', + fn=2; + otherwise, + fn=3; + %error(['unrecognized scoring fn ' scoring_fn]); + end + + if ~cache_full + place=cache(1,1); + else + [ignore place]=max(rand(1,L)); place=place+1; + end + + cache(place,:)=0; + cache(place,ps)=1; + cache(place,N+1)=j; + cache(place,N+2)=score; + cache(place,N+3)=fn; + + cache(1,1)=place+1; + if place>L || cache(1,2)~=0 + cache(1,2)=1; + end +end + + +%========================================================================================= +function [bool, score] = score_find_in_cache(cache,j,ps,scoring_fn) +% cache = score_find_in_cache(cache,j,ps,scoring_fn) +% +% V1.1 : 5 may 2003 (O. Francois, Ph. Leray) + + +%tic +L=size(cache,1)-1; +N=size(cache,2)-3; + +if N<1 + bool=0; + score=0; + return +end + +parents=zeros(1,N); +parents(ps)=1; +%parents(N+1)=j; + +switch scoring_fn +case 'bic', + fn=1; +case 'bayesian', + fn=2; +otherwise, + fn=3; + %error(['unrecognized scoring fn ' scoring_fn]); +end + +tmp=find(cache(2:L+1,N+3)==fn); +tmp=tmp+1; +tmp2=find(cache(tmp,N+1)==j); +candidats=tmp(tmp2); + +i=1; +while i<=N & ~isempty(candidats) + tmp=find(cache(candidats,i)==parents(i)); + candidats=candidats(tmp); + i=i+1; +end + +%Tpre=toc + +bool=~isempty(candidats); + +if bool + score=cache(candidats(1),N+2); +else + score=0; +end diff --git a/sourcecodes/bnt-master/SLP/scoring/score_find_in_cache.m b/sourcecodes/bnt-master/SLP/scoring/score_find_in_cache.m new file mode 100644 index 00000000..61f282e9 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/score_find_in_cache.m @@ -0,0 +1,63 @@ +function [bool, score] = score_find_in_cache(cache,j,ps,scoring_fn) +% [bool, score] = score_find_in_cache(cache,j,ps,scoring_fn) +% +% francois.olivier.c.h@gmail.com + +%tic +L=size(cache,1); +N=size(cache,2)-3; + +if N<1 + bool=0; + score=0; + return +end + +parents=zeros(1,N+1); +parents(ps)=1;parents(N+1)=j; + +switch scoring_fn + case 'bic', + fn=1; + case 'bayesian', + fn=2; + otherwise, + fn=3; + %error(['unrecognized scoring fn ' scoring_fn]); +end + +%parent = str2num(num2str(parents,'%1d')); +%[tmp y]=find(cache(:,N+3)==fn); +%if ~isempty(tmp) +% [tmp2 y]=find(str2num(num2str(cache(tmp,1:N+1),'%1d'))==parent); +% candidats=tmp(tmp2); +%else +% candidats=[]; +%end + +[tmp y]=find(cache(2:L,N+3)==fn); +tmp=tmp+1; +[tmp2 y]=find(cache(tmp,N+1)==j); +candidats=tmp(tmp2); +if ~isempty(candidats) + for i=1:N % N=size(cache,2)-3; + if ~isempty(candidats) + [tmp2 y]=find(cache(tmp,i)==parents(i)); + candidats=intersect(candidats,tmp(tmp2)); + end + end +end + +%Tpre=toc + +if ~isempty(candidats) + bool=1; +else + bool=0; +end + +if bool + score=cache(candidats(1),N+2); +else + score=0; +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/scoring/score_init_cache.m b/sourcecodes/bnt-master/SLP/scoring/score_init_cache.m new file mode 100644 index 00000000..fcbca268 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/scoring/score_init_cache.m @@ -0,0 +1,29 @@ +function cache = score_init_cache(N,L) +% SCORE_INIT_CACHE generate an empty cache for local computation in structure learning +% cache = score_init_cache(number_of_nodes,cache_size) +% +% For 2 nodes with cache of size 5 : +% +% cache = +% Nw b 0 0 0 --> Nw=number of writings in cache (+1) and b==1 iff the cache is full +% 0 0 1 -239.12 1 --> 1st familly in the cache (node 1 without parents) calculate with bic +% 0 0 2 -318.98 1 +% 1 0 2 -189.23 2 --> 3rd familly in the cache (node 2 with 1 as parent) calculate with bayesian +% 0 1 1 -251.09 1 +% 0 0 0 0 0 --> empty entry +% | | | | | +% | | | | |___> scoring function : 1 for 'bic', 2 for 'bayesian', ... +% | | | |__________> local score of the familly +% | | |_________________> son node of the familly +% | |__________________________> ==1 iff node 2 is parent of son node +% |______________________________> ==1 iff node 1 is parent of son node +% +% +% V1.1 : 6 may 2003 (O. Francois - francois.olivier.c.h@gmail.com, Ph. Leray - philippe.leray@univ-nantes.fr) +% +% + +cache=zeros(L+1,N+3); +cache(1,1)=2; + +% using a sparse matrix does not improve performances -- cgit 1.4.1