diff options
Diffstat (limited to 'sourcecodes/bnt-master/SLP/scoring')
19 files changed, 1733 insertions, 0 deletions
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 <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 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 <dhoiem@cs.cmu.edu>) + +[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 |