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