diff options
| author | ziejd2 | 2017-09-28 15:04:40 -0500 |
|---|---|---|
| committer | ziejd2 | 2017-09-28 15:04:40 -0500 |
| commit | 8070dc963753142bb86c4ed698d91fd623ed28e7 (patch) | |
| tree | d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/SLP/misc | |
| parent | 7cc31810d53176e805532b2789955f4eedbce6bb (diff) | |
| download | BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz | |
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
Diffstat (limited to 'sourcecodes/bnt-master/SLP/misc')
43 files changed, 4147 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/SLP/misc/CPT_from_bnet.m b/sourcecodes/bnt-master/SLP/misc/CPT_from_bnet.m new file mode 100644 index 00000000..64a62516 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/CPT_from_bnet.m @@ -0,0 +1,16 @@ +function CPT = CPT_from_bnet(bnet,v) +% Export CPT from bnet +% CPT = CPT_from_bnet(bnet,v) +% +% optional : v~=0 --> verbose mode + +if nargin<2, v=0; end +N = size(bnet.dag,1); +CPT = cell(1,N); +for j=1:N + CPD=struct(bnet.CPD{j}); + %counts{j}=CPD.counts; + CPT{j}=CPD.CPT; +end +if v, celldisp(CPT);end +%nsamples = CPD.nsamples; diff --git a/sourcecodes/bnt-master/SLP/misc/Markov_equivalent_dags.m b/sourcecodes/bnt-master/SLP/misc/Markov_equivalent_dags.m new file mode 100644 index 00000000..455037c9 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/Markov_equivalent_dags.m @@ -0,0 +1,42 @@ +function [n_dags, dag_list] = Markov_equivalent_dags(dag) + +% +% [n_dags, dag_list] = Markov_equivalent_dags(dag) +% +% generates a cell array of all Markov equivalent DAGs +% corresponding to the input DAG. +% +% YOU NEED TO HAVE THE STRUCTURE LEARNING PACKAGE IN PLACE TO USE THIS FUNCTION! +% +% Input: DAG (in standard format, i.e. dag(a,b)=1 if and only if a->b) +% +% Output: Number of DAGs generated and +% Cell array of all Markov-equivalent DAGs (in same format as input) +% +% Sample Use: +% +% Example 1: +% % Find all DAGs equivalent to DAG of Asia Network +% BN = mk_asia_bnet(); +% dag = BN.dag; +% [n_dags, dag_list] = Markov_equivalent_dags(dag); +% n_dags % Answer should be: 6 +% dag_list{1} % displays the first DAG, etc. +% +% Example 2: +% % Find all DAGs equivalent to random DAG +% dag = mk_rnd_dag(4); +% [n_dags, dag_list] = Markov_equivalent_dags(dag); +% +% Imme Ebert-Uphoff (ebert@tree.com), 2007 +% + + % find completed PDAG corresponding to DAG + cpdag = dag_to_cpdag(dag); + + % convert to our notation, i.e. directed edge has (-1) instead of (1) + signed_pdag = pdag_unsigned_to_signed(cpdag); + + % find all corresponding DAGs + [n_dags,dag_list] = pdag_to_all_dags( signed_pdag ); + diff --git a/sourcecodes/bnt-master/SLP/misc/bnt_to_mat.m b/sourcecodes/bnt-master/SLP/misc/bnt_to_mat.m new file mode 100644 index 00000000..27608feb --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/bnt_to_mat.m @@ -0,0 +1,18 @@ +function res = bnt_to_mat(data,misv) +% res = bnt_to_mat(data,misv) +% + +if nargin <1, error('Requires at least 1 argument.'), end +if nargin == 1, misv=-9999; end +taille=size(data); +long=taille(1); +larg=taille(2); +for i=1:long + for j=1:larg + if ~isempty(data{i,j}) + res(i,j)=data{i,j}; + else + res(i,j)=misv; + end + end +end diff --git a/sourcecodes/bnt-master/SLP/misc/chi2_table.m b/sourcecodes/bnt-master/SLP/misc/chi2_table.m new file mode 100644 index 00000000..e753d126 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/chi2_table.m @@ -0,0 +1,115 @@ +function value = chi2_table(P, df) +% function value = chi2_table(P, df) +% +% df 0.25 0.20 0.15 0.10 0.05 0.025 0.02 0.01 0.005 0.0025 0.001 0.0005 +% 1 1.32 1.64 2.07 2.71 3.84 5.02 5.41 6.63 7.88 9.14 10.83 12.12 +% 2 2.77 3.22 3.79 4.61 5.99 7.38 7.82 9.21 10.60 11.98 13.82 15.20 +% 3 4.11 4.64 5.32 6.25 7.81 9.35 9.84 11.34 12.84 14.32 16.27 17.73 +% 4 5.39 5.59 6.74 7.78 9.49 11.14 11.67 13.23 14.86 16.42 18.47 20.00 +% 5 6.63 7.29 8.12 9.24 11.07 12.83 13.33 15.09 16.75 18.39 20.51 22.11 +% 6 7.84 8.56 9.45 10.64 12.53 14.45 15.03 16.81 13.55 20.25 22.46 24.10 +% 7 9.04 5.80 10.75 12.02 14.07 16.01 16.62 18.48 20.28 22.04 24.32 26.02 +% 8 10.22 11.03 12.03 13.36 15.51 17.53 18.17 20.09 21.95 23.77 26.12 27.87 +% 9 11.39 12.24 13.29 14.68 16.92 19.02 19.63 21.67 23.59 25.46 27.83 29.67 +% 10 12.55 13.44 14.53 15.99 18.31 20.48 21.16 23.21 25.19 27.11 29.59 31.42 +% 11 13.70 14.63 15.77 17.29 19.68 21.92 22.62 24.72 26.76 28.73 31.26 33.14 +% 12 14.85 15.81 16.99 18.55 21.03 23.34 24.05 26.22 28.30 30.32 32.91 34.82 +% 13 15.93 15.58 18.90 19.81 22.36 24.74 25.47 27.69 29.82 31.88 34.53 36.48 +% 14 17.12 18.15 19.4 21.06 23.68 26.12 26.87 29.14 31.32 33.43 36.12 38.11 +% 15 18.25 19.31 20.60 22.31 25.00 27.49 28.26 30.58 32.80 34.95 37.70 39.72 +% 16 19.37 20.47 21.79 23.54 26.30 28.85 29.63 32.00 34.27 36.46 39.25 41.31 +% 17 20.49 21.61 22.98 24.77 27.59 30.19 31.00 33.41 35.72 37.95 40.79 42.88 +% 18 21.60 22.76 24.16 25.99 28.87 31.53 32.35 34.81 37.16 39.42 42.31 44.43 +% 19 22.72 23.90 25.33 27.20 30.14 32.85 33.69 36.19 38.58 40.88 43.82 45.97 +% 20 23.83 25.04 26.50 28.41 31.41 34.17 35.02 37.57 40.00 42.34 45.31 47.50 +% 21 24.93 26.17 27.66 29.62 39.67 35.48 36.34 38.93 41.40 43.78 46.80 49.01 +% 22 26.04 27.30 28.82 30.81 33.92 36.78 37.66 40.29 42.80 45.20 48.27 50.51 +% 23 27.14 28.43 29.98 32.01 35.17 38.08 38.97 41.64 44.18 46.62 49.73 52.00 +% 24 28.24 29.55 31.13 33.20 36.42 39.36 40.27 42.98 45.56 48.03 51.18 53.48 +% 25 29.34 30.68 32.28 34.38 37.65 40.65 41.57 44.31 46.93 49.44 52.62 54.95 +% 26 30.43 31.79 33.43 35.56 38.89 41.92 42.86 45.64 48.29 50.83 54.05 56.41 +% 27 31.53 32.91 34.57 36.74 40.11 43.19 44.14 46.96 49.64 52.22 55.48 57.86 +% 28 32.62 34.03 35.71 37.92 41.34 44.46 45.42 48.28 50.99 53.59 56.89 59.30 +% 29 33.71 35.14 36.85 39.09 42.56 45.72 46.69 49.59 52.34 54.97 58.30 60.73 +% 30 34.80 36.25 37.99 40.26 43.77 46.98 47.96 50.89 53.67 56.33 59.70 62.16 +% 40 45.62 47.27 49.24 51.81 55.76 59.34 60.44 63.69 66.77 69.70 73.40 76.09 +% 50 56.33 53.16 60.35 63.17 67.50 71.42 72.61 76.15 79.49 82.66 86.66 89.56 +% 60 66.98 68.97 71.34 74.40 79.08 83.30 84.58 88.38 91.95 95.34 99.61 102.7 +% 80 88.13 90.41 93.11 96.58 101.9 106.6 108.1 112.3 116.3 120.1 124.8 128.3 +% 100 109.1 111.7 114.7 118.5 124.3 129.6 131.1 135.8 140.2 144.3 149.4 153.2 +% df 0.25 0.20 0.15 0.10 0.05 0.025 0.02 0.01 0.005 0.0025 0.001 0.0005 +% + +if nargin<2, error('not enought arguments'); end +if P>0.25, error('first argument too big'); end +if P<0.0005, error('first argument too low'); end +if df>100, error('second argument too big'); end +if df<1, error('second argument too low'); end + +possible_P = [0.25,0.20,0.15,0.10,0.05,0.025,0.02 ,0.01,0.005,0.0025,0.001,0.0005]; +possible_df = [1:30, 40 50 60 80 100]; + +if ~ismember(P,possible_P), + res = sort([possible_P P]); res=res(end:-1:1); + [tmp indxP]=find(res==P); + indxP1 = indxP-1; indxP2=indxP; + P2=[possible_P(indxP1), possible_P(indxP2)]; + force_P = normalise(abs(P2-P)); +else + P2=[P P]; force_P=[1,0];[tmp,indxP1]=find(possible_P==P);indxP2=indxP1+1; +end +if ~ismember(df,possible_df), + res = sort([possible_df df]); + [tmp indxdf]=find(res==df); + indxdf1 = indxdf-1; indxdf2=indxdf; + df2=[possible_df(indxdf1), possible_df(indxdf2)]; + force_df = normalise(abs(df2-df)); +else + df2=[df df]; force_df=[1,0];[tmp,indxdf1]=find(possible_df==df);indxdf2=indxdf1+1; +end + +table = [[1.32,1.64,2.07,2.71,3.84,5.02,5.41,6.63 ,7.88,9.14,10.83,12.12];... +[2.77,3.22,3.79,4.61,5.99,7.38,7.82,9.21 ,10.60,11.98,13.82,15.20 ];... +[4.11,4.64,5.32,6.25,7.81,9.35,9.84,11.34 ,12.84,14.32,16.27,17.73 ];... +[5.39,5.59,6.74,7.78,9.49,11.14,11.67,13.23 ,14.86,16.42,18.47,20.00 ];... +[6.63,7.29,8.12,9.24,11.07,12.83,13.33,15.09 ,16.75,18.39,20.51,22.11 ];... +[7.84,8.56,9.45,10.64,12.53,14.45,15.03,16.81 ,13.55,20.25,22.46,24.10 ];... +[9.04,5.80,10.75,12.02,14.07,16.01,16.62 ,18.48,20.28,22.04,24.32,26.02 ];... +[10.22,11.03,12.03,13.36,15.51,17.53,18.17 ,20.09,21.95,23.77,26.12,27.87 ];... +[11.39,12.24,13.29,14.68,16.92,19.02,19.63 ,21.67,23.59,25.46,27.83,29.67 ];... +[12.55,13.44,14.53,15.99,18.31,20.48,21.16 ,23.21,25.19,27.11,29.59,31.42 ];... +[13.70,14.63,15.77,17.29,19.68,21.92,22.62,24.72,26.76,28.73,31.26,33.14 ];... +[14.85,15.81,16.99,18.55,21.03,23.34,24.05,26.22,28.30,30.32,32.91,34.82 ];... +[15.93,15.58,18.90,19.81,22.36,24.74,25.47,27.69,29.82,31.88,34.53,36.48 ];... +[17.12,18.15,19.4,21.06,23.68,26.12,26.87,29.14,31.32,33.43,36.12,38.11 ];... +[18.25,19.31,20.60,22.31,25.00,27.49,28.26,30.58,32.80,34.95,37.70,39.72 ];... +[19.37,20.47,21.79,23.54,26.30,28.85,29.63,32.00,34.27,36.46,39.25,41.31 ];... +[20.49,21.61,22.98,24.77,27.59,30.19,31.00,33.41,35.72,37.95,40.79,42.88 ];... +[21.60,22.76,24.16,25.99,28.87,31.53,32.35,34.81,37.16,39.42,42.31,44.43 ];... +[22.72,23.90,25.33,27.20,30.14,32.85,33.69,36.19,38.58,40.88,43.82,45.97 ];... +[23.83,25.04,26.50,28.41,31.41,34.17,35.02,37.57,40.00,42.34,45.31,47.50 ];... +[24.93,26.17,27.66,29.62,39.67,35.48,36.34,38.93,41.40,43.78,46.80,49.01 ];... +[26.04,27.30,28.82,30.81,33.92,36.78,37.66,40.29,42.80,45.20,48.27,50.51 ];... +[27.14,28.43,29.98,32.01,35.17,38.08,38.97,41.64,44.18,46.62,49.73,52.00 ];... +[28.24,29.55,31.13,33.20,36.42,39.36,40.27,42.98,45.56,48.03,51.18,53.48 ];... +[29.34,30.68,32.28,34.38,37.65,40.65,41.57,44.31,46.93,49.44,52.62,54.95 ];... +[30.43,31.79,33.43,35.56,38.89,41.92,42.86,45.64,48.29,50.83,54.05,56.41 ];... +[31.53,32.91,34.57,36.74,40.11,43.19,44.14,46.96,49.64,52.22,55.48,57.86 ];... +[32.62,34.03,35.71,37.92,41.34,44.46,45.42,48.28,50.99,53.59,56.89,59.30 ];... +[33.71,35.14,36.85,39.09,42.56,45.72,46.69,49.59,52.34,54.97,58.30,60.73 ];... +[34.80,36.25,37.99,40.26,43.77,46.98,47.96,50.89,53.67,56.33,59.70,62.16 ];... +[45.62,47.27,49.24,51.81,55.76,59.34,60.44,63.69,66.77,69.70,73.40,76.09 ];... +[56.33,53.16,60.35,63.17,67.50,71.42,72.61,76.15,79.49,82.66,86.66,89.56 ];... +[66.98,68.97,71.34,74.40,79.08,83.30,84.58,88.38,91.95,95.34,99.61,102.7 ];... +[88.13,90.41,93.11,96.58,101.9,106.6,108.1,112.3,116.3,120.1,124.8,128.3 ];... +[109.1,111.7,114.7,118.5,124.3,129.6,131.1,135.8,140.2,144.3,149.4,153.2]]; + +values = table(indxdf1:indxdf2,indxP1:indxP2); +value1 = values(1,:)*force_P'; +value2 = values(2,:)*force_P'; +value3 = values(:,1)'*force_df'; +value4 = values(:,2)'*force_df'; +value5 = [value1 value2]*force_df'; +value6 = [value3 value4]*force_P'; +value = mean([value5, value6]); + diff --git a/sourcecodes/bnt-master/SLP/misc/classification_evaluation.m b/sourcecodes/bnt-master/SLP/misc/classification_evaluation.m new file mode 100644 index 00000000..a2b107bc --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/classification_evaluation.m @@ -0,0 +1,35 @@ +function [ratio, ratiominus, ratioplus, proba_post, yt] = classification_evaluation(bnet, BDT, class) +% Computes the classification ratio of a bnet structure on a test dataset BDT +% [ratio ratiominus ratioplus] = classification_evaluation(bnet, BDT, class) +% +% [ratiominus rationplus] is the 95 percent confident interval. +% results are in percentage [0 100]. +% +% francois.olivier.c.h@gmail.com +% + + [proba_post,engine] = inference(bnet, mat_to_bnt(BDT), class); + [tmp yt] = max(proba_post, [],2); + [N L] = size(BDT); + count = length(find(BDT(class,:)==yt')); + ratio = 100*count/L; +[ratiominus, ratioplus] = confiance(ratio,L); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [I, J] = confiance(t, N) +% Compute the 95 percent confident interval +% +% see Y. Bennani and F. Bossaert, +% Predictive neural networks for traffic disturbance detection in the telephone network +% In Proceedings of IMACS-CESA 1996, Lille, France. + +Z = 1.96; % this value for the 95 percent confident interval +T = t/100; +tmp = (Z*Z)/N; +D = 1+tmp; +N1 = T+tmp/2; +tmp2 = T*(1-T)/N + tmp/(4*N); +N2 = Z*sqrt(tmp2); +I = 100*(N1-N2)/D; +J = 100*(N1+N2)/D; diff --git a/sourcecodes/bnt-master/SLP/misc/complete_pattern.m b/sourcecodes/bnt-master/SLP/misc/complete_pattern.m new file mode 100644 index 00000000..74cb8db7 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/complete_pattern.m @@ -0,0 +1,119 @@ +function completed_pdag = complete_pattern(pdag) + +% +% completed_dag = complete_pattern(pdag) +% +% uses Rules R1-R4 of Meek (1995) to complete +% orientations in a pdag as far as possible, +% i.e. every compelled edge is oriented. +% +% (Rules R1-R4 are also summarized in Pearl (2000), p.51 +% and Neapolitan (2004), p. 546.) +% +% Since the PC algorithm also uses Rules R1-R3, their implementation was +% copied (with some modifications) from function learn_struct_pdag_pc. +% +% Rule R4 is necessary here, since the orientations in the input +% pdag do not just represent v-structures. +% +% Imme Ebert-Uphoff (ebert@tree.com), 2007 +% + + DIAGNOSTICS_ON = false; + + n = length(pdag); + old_pdag = zeros(n); + %iter = 0; + while ~isequal(pdag, old_pdag) + %iter = iter + 1; + old_pdag = pdag; + + % Rule R1 + [A,B] = find(pdag==-1); % a -> b + for i=1:length(A) + a = A(i); b = B(i); + undirected = abs(pdag) + abs(pdag)'; + % Adjacency test in undirected matrix: + % a adjacent b <=> undirected(a,b) ==0 + % That's easier to use than adjacency test in pdag: + % a adjacent b <=> pdag(a,b)==0 and pdag(b,a)==0 + + % Find all nodes c such that b-c and c not adjacent a + C = find(pdag(b,:)==1 & undirected(a,:)==0); + if ~isempty(C) + pdag(b,C) = -1; pdag(C,b) = 0; + if DIAGNOSTICS_ON + for j=1:length(C) + fprintf('Rule 1: %d -> %d\n', b, C(j)); + end + end + end + end + + % Rule R2 + [A,B] = find(pdag==1); % unoriented a-b edge + for i=1:length(A) + a = A(i); b = B(i); + if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' ); + pdag(a,b) = -1; pdag(b,a) = 0; + if DIAGNOSTICS_ON + fprintf('Rule 2: %d -> %d\n', a, b); + end + end + end + + % Rule R3 + [A,B] = find(pdag==1); % a-b + for i=1:length(A) + a = A(i); b = B(i); + C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' ); + % C contains nodes c s.t. a-c->b-a + + % Extract lines and columns corresponding only to the set of nodes C + core = pdag(C,C); + + % Prepare adjacency test: + unoriented = abs(core) + abs(core)'; + % Now: a non-adjacent b <==> unoriented(a,b) == 0 + + % Prepare to detect existence of non-adjacent pairs of nodes in C. + % Set diagonal to 1, to prevent finding pairs of IDENTICAL nodes: + unoriented = setdiag(unoriented, 1); + if any(unoriented(:)==0) % C contains 2 different non adjacent elements + pdag(a,b) = -1; pdag(b,a) = 0; + if DIAGNOSTICS_ON + fprintf('Rule 3: %d -> %d\n', a, b); + end + end + end + + % Rule 4 + [A,B] = find(pdag==1); % unoriented a-b edge + for i=1:length(A) + a = A(i); b = B(i); + + % Prepare adjacency test: + % unoriented(i,j) is 0 (non-adj) or 1 (directed) or 2 (undirected) + unoriented = abs(pdag) + abs(pdag)'; + + % Find c such that c -> b and a,c are adjacent (a-c or a->c or a<-c) + C = find( (pdag(:,b)==-1)' & (unoriented(a,:)>=1) ); + for j=1:length(C) + c = C(j); + % Check whether there is any node d, such that + % d->c AND a-d AND b NOT adjacent to d + if any( (pdag(:,c)==-1)' & (pdag(a,:)==1) & (unoriented(b,:)==0) ) + pdag(a,b) = -1; pdag(b,a) = 0; + if DIAGNOSTICS_ON + fprintf('Rule 4: %d -> %d\n', a, b); + end + end + end + end + + end % end of while + + % Oriented all possible edges. Return result. + completed_pdag = pdag; +end + diff --git a/sourcecodes/bnt-master/SLP/misc/compute_bnet_nparams.m b/sourcecodes/bnt-master/SLP/misc/compute_bnet_nparams.m new file mode 100644 index 00000000..3366375f --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/compute_bnet_nparams.m @@ -0,0 +1,13 @@ +function [D,d]= compute_bnet_nparams(bnet) +% [D,d] = compute_bnet_nparams(bnet) +% +% D is the dimension of the network +% d is the vector containing the number of parameters for all nodes + +N = length(bnet.dag); +d = zeros(1,N); +for i=1:N + a = struct(bnet.CPD{i}); + d(i) = a.nparams; +end +D = sum(d); diff --git a/sourcecodes/bnt-master/SLP/misc/confiance.m b/sourcecodes/bnt-master/SLP/misc/confiance.m new file mode 100644 index 00000000..a9416bf8 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/confiance.m @@ -0,0 +1,16 @@ +function [I, J] = confiance(t, N) +% Compute the 95 percent confident interval +% +% see Y. Bennani and F. Bossaert, +% Predictive neural networks for traffic disturbance detection in the telephone network +% In Proceedings of IMACS-CESA 1996, Lille, France. + +Z = 1.96; % this value for the 95 percent confident interval +T = t/100; +tmp = (Z*Z)/N; +D = 1+tmp; +N1 = T+tmp/2; +tmp2 = T*(1-T)/N + tmp/(4*N); +N2 = Z*sqrt(tmp2); +I = 100*(N1-N2)/D; +J = 100*(N1+N2)/D; diff --git a/sourcecodes/bnt-master/SLP/misc/cpdag_to_dag.m b/sourcecodes/bnt-master/SLP/misc/cpdag_to_dag.m new file mode 100644 index 00000000..8e106610 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/cpdag_to_dag.m @@ -0,0 +1,52 @@ +function dag = cpdag_to_dag(cpdags) +% dags = cpdag_to_dag(cpdags) +% +% CPDAG_TO_DAG produce a N*N matrix of a dag which instantiate cpdag. +% (also works with a cell array of cpdags, returning a cell array of dags) +% make sur that your entry is a completed PDAG +% this function can't be use instead of PDAG_TO_DAG +% +% see Chickering (2002) : Learning equivalence classes of bayesian networks, JMLR2, pp475-479 +% +% francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr +% 31 march 2006 + +if ~iscell(cpdags) + cpdag=cell(1,1); + cpdag{1}=cpdags; +else + cpdag=cpdags; +end + +for da=1:length(cpdag) + + N=length(cpdag{da}); + dag=cpdag{da}; + + unprocessed = find_nodes_in_undirected_component(dag); + while ~isempty(unprocessed) + nbr_parents = []; + for i=1:length(unprocessed) + nbr_parents(end+1)=length(parents(dag-dag.*dag',unprocessed(i))); %nbr_parents(end+1)=length(parents(dag,unprocessed(i))); + end + [tmp, idx] = max(nbr_parents); + node = unprocessed(idx); + dag(parents(dag.*dag',node),node)=0; %dag(parents(dag,node),node)=0; + %dag(myintersect(parents(dag,node),unprocessed), node)=0; % Wei Lu + unprocessed = find_nodes_in_undirected_component(dag); + end + + dags{da}=dag; +end + +if ~iscell(cpdags) + dag=dags{1}; +else + dag=dags; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function unprocessed = find_nodes_in_undirected_component(dag) +undirected_edges = dag.*dag'; +[unprocessed, tmp] = find(undirected_edges); +unprocessed = unique(unprocessed); diff --git a/sourcecodes/bnt-master/SLP/misc/cpdag_to_dag2.m b/sourcecodes/bnt-master/SLP/misc/cpdag_to_dag2.m new file mode 100644 index 00000000..1057a9eb --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/cpdag_to_dag2.m @@ -0,0 +1,52 @@ +function dag = cpdag_to_dag2(cpdags) +% dags = cpdag_to_dag(cpdags) +% +% CPDAG_TO_DAG produce a N*N matrix of a dag which instantiate cpdag. +% (also works with a cell array of cpdags, returning a cell array of dags) +% make sur that your entry is a completed PDAG +% this function can't be use instead of PDAG_TO_DAG +% +% francois.olivier.c.h@gmail.com +% 7 may 2003 - OLD version + +if ~iscell(cpdags) + cpdag=cell(1,1); + cpdag{1}=cpdags; +else + cpdag=cpdags; +end + +for da=1:length(cpdag) + + N=length(cpdag{da}); + dag=cpdag{da}; dag2=dag; + unprocessed = []; + + for i=1:(N-1) + for j=(i+1):N + if dag2(i,j)==1 & dag2(j,i)==1 + if ~myismember(i,unprocessed) + unprocessed = [unprocessed;i]; + end + if ~myismember(j,unprocessed) + unprocessed = [unprocessed;j]; + end + end + end + end + + for i=1:length(unprocessed) + for j=1:N + if dag(unprocessed(i),j)==1 & dag(j,unprocessed(i))==1 + dag(j,unprocessed(i))=0; + end + end + end + dags{da}=dag; +end + +if ~iscell(cpdags) + dag=dags{1}; +else + dag=dags; +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/dag_to_cpdag.m b/sourcecodes/bnt-master/SLP/misc/dag_to_cpdag.m new file mode 100644 index 00000000..70ec3a7e --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/dag_to_cpdag.m @@ -0,0 +1,115 @@ +function [cpdag] = dag_to_cpdag(dags) +% (also works with a cell array of dags, returning a cell array of cpdags) +% DAG_TO_CPDAG produce a N*N matrix which values respect : +% +% If the edge is compelled then 1 on the edge. +% If the edge is reversible then 1 on the edge and 1 in the reverse edge. +% +% Make sure that the entry is a DAG. +% +% See D.M. Chickering: "Learning Equivalence Classes of Bayesian-Network Structures". +% +% +% francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr, alain.delaplace@univ-tours.fr + +if ~iscell(dags) + dag=cell(1,1); + dag{1}=dags; +else + dag=dags; +end + +for da=1:length(dag) + cpdags{da} = abs(label_edges(dag{da})); +end + +if ~iscell(dags) + cpdag=cpdags{1}; +else + cpdag=cpdags; +end + +%%============================================================================== + +function [label] = label_edges(dag) +% LABEL-EDGES produce a N*N matrix which values are +% +1 if the edge is compelled or +% -1 if the edge is reversible. +% Make sure that the entry is a DAG. +% +% francois.olivier.c.h@gmail.com + +N=length(dag); +[order xedge yedge] = order_edges(dag); +label = 2*dag; + +NbEdges = length(xedge) ; + +for Edge=1:NbEdges, + xlow=xedge(Edge); + ylow=yedge(Edge); + if label(xlow,ylow)==2 + fin = 0; + wcompelled = find(label(:,xlow)==1); + parenty = find(label(:,ylow)~=0); + + %for w = wcompelled + for s = 1:length(wcompelled) + w = wcompelled(s); + if ~ismember(w,parenty) + label(parenty,ylow)=1; + label(ylow,parenty)=0; + fin = 1; + elseif fin == 0 + label(w,ylow)=1; + label(ylow,w)=0; % + end + end + if fin == 0 + parentx = [xlow ; find(label(:,xlow)~=0)]; + if ~isempty(mysetdiff(parenty,parentx)) + label(xlow,ylow)=1; % + label(ylow,xlow)=0; % + ttp=find(label(:,ylow)==2); + label(ttp,ylow)=1; + label(ylow,ttp)=0; % + else + ttp=find(label(:,ylow)==2); + label(ttp,ylow)=-1; + label(ylow,ttp)=-1; + end + end + end +end + +%%%======================================================================================== +function [order, x, y] = order_edges(dag) +% ORDER_EDGES produce a total (natural) ordering over the edges in a DAG. +% Make sure that the entry is a DAG. +% +% francois.olivier.c.h@gmail.com +% +% 2 mai 2003 + +if acyclic(dag)==0 + error('Requires an acyclic graph'); +end + +N=length(dag); +order = zeros(N,N); + +node_order = topological_sort(dag); +[tmp oo] = sort(node_order); + +dag=dag(node_order,node_order); +[x y]=find(flipud(dag)==1); +nb_edges=length(x); + +if nb_edges~=0 + order(sub2ind([N N],N+1-x,y))=1:nb_edges ; +end + +order=order(oo,oo); +x=node_order(N+1-x); +y=node_order(y); + diff --git a/sourcecodes/bnt-master/SLP/misc/dag_to_cpdag1.m b/sourcecodes/bnt-master/SLP/misc/dag_to_cpdag1.m new file mode 100644 index 00000000..c3189a1d --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/dag_to_cpdag1.m @@ -0,0 +1,128 @@ +function [cpdag] = dag_to_cpdag1(dags) +% 2 +% (also works with a cell array of dags, returning a cell array of cpdags) +% DAG_TO_CPDAG produce a N*N matrix which values respect : +% +% If the edge is compelled then 1 on the edge. +% If the edge is reversible then 1 on the edge and 1 in the reverse edge. +% +% Make sure that the entry is a DAG. +% +% See D.M. Chickering: "Learning Equivalence Classes of Bayesian-Network Structures". +% +% +% francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr + +if ~iscell(dags) + dag=cell(1,1); + dag{1}=dags; +else + dag=dags; +end + +for da=1:length(dag) + cpdags{da} = abs(label_edges(dag{da})); +end + +if ~iscell(dags) + cpdag=cpdags{1}; +else + cpdag=cpdags; +end + +%%============================================================================== + +function [label] = label_edges(dag) +% LABEL-EDGES produce a N*N matrix which values are +% +1 if the edge is compelled or +% -1 if the edge is reversible. +% Make sure that the entry is a DAG. +% +% francois.olivier.c.h@gmail.com + +N=length(dag); +[order xedge yedge] = order_edges(dag); +label = 2*dag; % all edges as unknown + +NbEdges = length(xedge) ; +%xedge=x, yedge=y, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +for Edge=1:NbEdges, + xlow=xedge(Edge); + ylow=yedge(Edge); + if label(xlow,ylow)==2 + fin = 0; + + %wcompelled = find(label(:,xlow)~=2); + wcompelled = find(label(:,xlow)==1); + + parenty = find(label(:,ylow)~=0); + %sonsy = find(label(ylow,:)~=0); + + for s = 1:length(wcompelled) + w = wcompelled(s); + if ~ismember(w,parenty) + label(xlow,ylow)=1; + label(ylow,xlow)=0; % + label(parenty,ylow)=1; + label(ylow,parenty)=0; % + %label(ylow,sonsy)=1; % + %label(sonsy,ylow)=0; % + fin = 1; + + elseif fin == 0 + label(w,ylow)=1; + %label(ylow,w)=0; % + end + end + if fin == 0 + parentx = [xlow ; find(label(:,xlow)~=0)]; + if ~isempty(mysetdiff(parenty,parentx)) + %label(xlow,ylow)=1; % + %label(ylow,xlow)=0; % + + label(find(label(:,ylow)==2),ylow)=1; + label(ylow,find(label(ylow,:)==2))=1; % + else + label(xlow,ylow)=-1; + label(ylow,xlow)=-1; % + ttp = find(label(:,ylow)==2); + label(ttp,ylow)=-1; + label(ylow,ttp)=-1; % + end + end + end +end + +%%%======================================================================================== +function [order, x, y] = order_edges(dag) +% ORDER_EDGES produce a total (natural) ordering over the edges in a DAG. +% Make sure that the entry is a DAG. +% +% francois.olivier.c.h@gmail.com +% +% 2 mai 2003 + +if acyclic(dag)==0 + error('Requires an acyclic graph'); +end + +N=length(dag); +order = zeros(N,N); + +node_order = topological_sort(dag); +[tmp oo] = sort(node_order); + +dag2=dag(node_order,node_order); +[x y]=find(flipud(dag2)==1); +nb_edges=length(x); + +if nb_edges~=0 + order(sub2ind([N N],N+1-x,y))=1:nb_edges ; +end + +order=order(oo,oo); +x=node_order(N+1-x); +y=node_order(y); + diff --git a/sourcecodes/bnt-master/SLP/misc/discretization.m b/sourcecodes/bnt-master/SLP/misc/discretization.m new file mode 100644 index 00000000..870583af --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/discretization.m @@ -0,0 +1,59 @@ +function [appD, testD, bornes] = discretization(critere, continues, miss, app, test) +% [appD, testD, bornes] = discretization(critere, continious, miss, app, test) +% +% Inputs : +% critere = 1, 2, 3 or 4 (see hict_ic for details) +% continious = vector of continious variables to discretize +% miss ~= 0 if it exists missing values coded by 'miss' value +% app = Learning base +% test = Test base (only the learning base is used to make the discretization rules) [optionnal] +% +% Outputs : +% appD = Learning base with discretized entries on 'continious' variables +% testD = Test base with discretized entries on 'continious' variables +% bornes = limits of discretization intervals found by hist_ic +% + +tt=1; +if nargin<5, test=[]; tt=0; end +app = app'; +test = test'; +[ma, Na] = size(app); +[mt, Nt] = size(test); +I = []; +testD = []; + +if miss, + [I J]=find(app==miss); + [I2 J2]=find(test==miss); +end +completes=setdiff(1:ma,I); +donnees_continue=app(completes,continues); + +% echantillonnage +[n,bornes,nbbornes,xx]=hist_ic(donnees_continue,critere); + +% on re-distribue l'ensemble des donnees d'apprentissage continues +[n2,appD_continues]=histc_ic(app(:,continues),bornes); +if tt, [n2test,testD_continues]=histc_ic(test(:,continues),bornes); end + +% on insere les donnees continues discretisees dans les matrices +appD=app; +if tt, testD=test; end +for i=1:length(continues) + appD(:,continues(i))=appD_continues(:,i); + if tt, testD(:,continues(i))=testD_continues(:,i); end +end + +if miss, + for k=1:length(I) + app(I(k),J(k))=miss; + end + if tt, + for l=1:length(I2) + testD(I2(l),J2(l))=miss; + end + end +end +appD = appD'; +testD = testD'; \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/editing_dist.m b/sourcecodes/bnt-master/SLP/misc/editing_dist.m new file mode 100644 index 00000000..5b697752 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/editing_dist.m @@ -0,0 +1,17 @@ +function D = editing_dist(dag1, dag2) +% D = editing_dist(dag1, dag2) +% +% d = 1 if arc1 <> arc2 +% + +[n1 m1]=size(dag1); +[n2 m2]=size(dag2); + +if n1~=m1 | n1~=n2 | n2~=m2 | m1~=m2 + error('formats non compatibles'); + D=-inf; +end + +de=abs(dag1-dag2); +mauv=find(triu(or(de,de'))); +D=size(mauv,1); diff --git a/sourcecodes/bnt-master/SLP/misc/export_cases.m b/sourcecodes/bnt-master/SLP/misc/export_cases.m new file mode 100644 index 00000000..4ab0d790 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/export_cases.m @@ -0,0 +1,57 @@ +function S = export_cases(data, names, file, misv) +% filepath = export_cases(data, names, 'filename', missing_value) +% filename and missing_value [default [] if iscell(data) or -9999 if not] are optional +% +% Exports BNT datasets to Netica cases (* for missing data) +% +% Written by Francois.Olivier.C.H@gmail.com +% +% Informations could be found here http://www.norsys.com/downloads/ +% +% version 80807 + +% inits +if nargin==1, file=['dnet' datestr(now,'-yymmdd-HHMMSS')]; end +if nargin<2, error('Variable Names needed'); end +if length(file)>5, + if prod(double(file((end-3):end)~='.cas')), file=[file '.cas']; end + name = file(1:end-5); +else + name = file; + file = [file '.cas']; +end +if nargin<4, misv=-9999; end +if iscell(data), data = bnt_to_mat(data,misv); end + +[N m]=size(data); +if length(names)~=N, error('Sizes must be the same'); end + +% header of the file +fid = fopen(file, 'w'); +fprintf(fid, '// exported from the Bayes Net Toolbox with export_cases function \n'); +fprintf(fid, '// please report bugs to francois.olivier.c.h@gmail.com\n\n'); + +% write names +for i=1:N + fprintf(fid, '%s\t',names{i}); +end + +% exports +for l=1:m, + fprintf(fid, '\n'); + for i=1:N, + if data(i,l)~=misv, + fprintf(fid, '%s',['x' num2str(data(i,l))]); + else + fprintf(fid, '%s','*'); + end + if i<N, fprintf(fid, '\t'); end + end +end + +% closes file +fprintf(fid,'\n'); +fclose(fid); + +% outputs string +S = [pwd '/' file]; diff --git a/sourcecodes/bnt-master/SLP/misc/export_dnet.m b/sourcecodes/bnt-master/SLP/misc/export_dnet.m new file mode 100644 index 00000000..9b0cbd04 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/export_dnet.m @@ -0,0 +1,153 @@ +function S = export_dnet(bnet, file, proba) +% filepath = export_dnet(bnet, 'filename', includeparameters) +% filename and includeparameters ([0] or 1) are optional +% +% Exports BNT bnets to Netica dnet +% +% Written by Francois.Olivier.C.H@gmail.com +% +% Informations could be found here http://www.norsys.com/downloads/ +% +% Version 80806 +% Supports only discret bayesian network + +if nargin==1, file=['dnet' datestr(now,'-yymmdd-HHMMSS')]; end +if nargin<3, proba=0; end + +% generating filename +if length(file)>4, + if prod(double(file((end-4):end)~='.dnet')), file=[file '.dnet']; end + name = file(1:end-5); +else + name = file; + file = [file '.dnet']; +end + +% generating node names 1:N if non existant +N=length(bnet.dag); +if isempty(bnet.names), + for i=1:N, + keys{i}=['' num2str(i) '']; + vals{i}=i; + end + bnet.names = assocarray(keys, vals); +end + +% header of the file +fid = fopen(file, 'w'); +fprintf(fid, '// ~->[DNET-1]->~\n\n'); +fprintf(fid, '// exported from the Bayes Net Toolbox with export_dnet function \n'); +fprintf(fid, '// please report bugs to francois.olivier.c.h@gmail.com\n'); +if proba, fprintf(fid, '// Take care ! Parents'' order isn''t the same in probability table comments\n'); end +fprintf(fid, ['\nbnet ' name ' {']); + +% main loop +for node = 1:N + name = get_key(bnet.names,node); + + fprintf(fid, ['\nnode ' name ' {']); + + fprintf(fid, '\n\tkind = NATURE;'); + + fprintf(fid, '\n\tdiscrete = TRUE;'); + +% states names are x1,x2:xsize(node) + fprintf(fid, '\n\tstates = ('); + for l=1:bnet.node_sizes(node), + fprintf(fid, 'x%d',l); + if l~=bnet.node_sizes(node), fprintf(fid, ', '); end + end, fprintf(fid, ');'); + +% declare parent in counter order to be coherent with prob section + fprintf(fid, '\n\tparents = ('); + par = bnet.parents{node}; + for l=length(par):-1:1 + fprintf(fid, '%s',get_key(bnet.names,par(l))); + if l~=1, fprintf(fid, ', '); end + end, fprintf(fid, ');'); + +% fill probs if proba==1 + if proba, + fprintf(fid, '\n\tprobs =\n'); + % inits + fprintf(fid, '\t//\t'); + for l = 1:bnet.node_sizes(node), + fprintf(fid, 'x%d\t', l); + end + fprintf(fid, '\t//'); + if ~isempty(par), + for l = unique(par), %right order this time because of the way both BNT and netica work + fprintf(fid, '\t%s', get_key(bnet.names,l)); + end + end + + % opens tab + fprintf(fid, '\n\t'); + for l = unique(par), + fprintf(fid, '('); + end + + % fullfils probs + CPT = CPT_from_bnet(bnet); + CPT = CPT{node}; + CPT=CPT(:); % good order whatever the node size ???? + if isempty(par), + fprintf(fid, '\t'); + for i=1:length(CPT)-1, + fprintf(fid, '%1.4f, ', CPT(i)); + end + fprintf(fid, '%1.4f);',CPT(end)); + else % if there are parents + endi=0; + parsiz = prod(bnet.node_sizes([par])); + parentstates = ones(1,length(par)); parentstates(end)=0; + for i=1:parsiz + % prints probas + fprintf(fid, '(\t'); + for j = 1:bnet.node_sizes(node) + prob = CPT(i+parsiz*(j-1)); + if j~=bnet.node_sizes(node), fprintf(fid, '%1.4f, ',prob); + else fprintf(fid, '%1.4f',prob); end + end + + % closes parenthesis if needed + endii=endi; + if i~=parsiz, + while endi>0, fprintf(fid, ')');endi=endi-1; end + fprintf(fid, '),\t//\t'); + else % close tab + for l = unique(par), fprintf(fid, ')'); end + fprintf(fid, ');\t//\t'); + end + + % prints node states + res=i; + for l=1:length(par)-1 + if mod(i+1, prod(bnet.node_sizes(par(1:l))))==0, endi=endi+1; end % counts parenthesis + + resaff = mod(res, bnet.node_sizes(l)); + if resaff==0, resaff=bnet.node_sizes(par(l)); end + res = div(res-1, bnet.node_sizes(par(l)))+1; + + fprintf(fid, '%d\t',resaff); + end + + resaff = div(i-1, prod(bnet.node_sizes(par(1:end-1))))+1; + fprintf(fid, '%d',resaff); + + % opens parenthesis if needed + if i~=parsiz, fprintf(fid, '\n'); end + fprintf(fid, '\t'); + if i~=parsiz, while endii>0, fprintf(fid, '(');endii=endii-1; end, end + end + end + end + fprintf(fid, '\n\t};'); +end + +% closes file +fprintf(fid,'\n};\n'); +fclose(fid); + +% outputs string +S = [pwd '/' file]; diff --git a/sourcecodes/bnt-master/SLP/misc/exportfig.m b/sourcecodes/bnt-master/SLP/misc/exportfig.m new file mode 100644 index 00000000..ef734fa3 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/exportfig.m @@ -0,0 +1,991 @@ +function varargout = exportfig(varargin) +%EXPORTFIG Export a figure. +% EXPORTFIG(H, FILENAME) writes the figure H to FILENAME. H is +% a figure handle and FILENAME is a string that specifies the +% name of the output file. +% +% EXPORTFIG(H, FILENAME, OPTIONS) writes the figure H to FILENAME +% with options initially specified by the structure OPTIONS. The +% field names of OPTIONS must be legal parameters listed below +% and the field values must be legal values for the corresponding +% parameter. Default options can be set in releases prior to R12 +% by storing the OPTIONS structure in the root object's appdata +% with the command +% setappdata(0,'exportfigdefaults', OPTIONS) +% and for releases after R12 by setting the preference with the +% command +% setpref('exportfig', 'defaults', OPTIONS) +% +% EXPORTFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies +% parameters that control various characteristics of the output +% file. Any parameter value can be the string 'auto' which means +% the parameter uses the default factory behavior, overriding +% any other default for the parameter. +% +% Format Paramter: +% 'Format' a string +% specifies the output format. Defaults to 'eps'. For a +% list of export formats type 'help print'. +% 'Preview' one of the strings 'none', 'tiff' +% specifies a preview for EPS files. Defaults to 'none'. +% +% Size Parameters: +% 'Width' a positive scalar +% specifies the width in the figure's PaperUnits +% 'Height' a positive scalar +% specifies the height in the figure's PaperUnits +% 'Bounds' one of the strings 'tight', 'loose' +% specifies a tight or loose bounding box. Defaults to 'tight'. +% 'Reference' an axes handle or a string +% specifies that the width and height parameters +% are relative to the given axes. If a string is +% specified then it must evaluate to an axes handle. +% +% Specifying only one dimension sets the other dimension +% so that the exported aspect ratio is the same as the +% figure's or reference axes' current aspect ratio. +% If neither dimension is specified the size defaults to +% the width and height from the figure's or reference +% axes' size. Tight bounding boxes are only computed for +% 2-D views and in that case the computed bounds enclose all +% text objects. +% +% Rendering Parameters: +% 'Color' one of the strings 'bw', 'gray', 'cmyk' +% 'bw' specifies that lines and text are exported in +% black and all other objects in grayscale +% 'gray' specifies that all objects are exported in grayscale +% 'rgb' specifies that all objects are exported in color +% using the RGB color space +% 'cmyk' specifies that all objects are exported in color +% using the CMYK color space +% 'Renderer' one of 'painters', 'zbuffer', 'opengl' +% specifies the renderer to use +% 'Resolution' a positive scalar +% specifies the resolution in dots-per-inch. +% 'LockAxes' one of 0 or 1 +% specifies that all axes limits and ticks should be fixed +% while exporting. +% +% The default color setting is 'bw'. +% +% Font Parameters: +% 'FontMode' one of the strings 'scaled', 'fixed' +% 'FontSize' a positive scalar +% in 'scaled' mode multiplies with the font size of each +% text object to obtain the exported font size +% in 'fixed' mode specifies the font size of all text +% objects in points +% 'DefaultFixedFontSize' a positive scalar +% in 'fixed' mode specified the default font size in +% points +% 'FontSizeMin' a positive scalar +% specifies the minimum font size allowed after scaling +% 'FontSizeMax' a positive scalar +% specifies the maximum font size allowed after scaling +% 'FontEncoding' one of the strings 'latin1', 'adobe' +% specifies the character encoding of the font +% 'SeparateText' one of 0 or 1 +% specifies that the text objects are stored in separate +% file as EPS with the base filename having '_t' appended. +% +% If FontMode is 'scaled' but FontSize is not specified then a +% scaling factor is computed from the ratio of the size of the +% exported figure to the size of the actual figure. +% +% The default 'FontMode' setting is 'scaled'. +% +% Line Width Parameters: +% 'LineMode' one of the strings 'scaled', 'fixed' +% 'LineWidth' a positive scalar +% 'DefaultFixedLineWidth' a positive scalar +% 'LineWidthMin' a positive scalar +% specifies the minimum line width allowed after scaling +% 'LineWidthMax' a positive scalar +% specifies the maximum line width allowed after scaling +% The semantics of 'Line' parameters are exactly the +% same as the corresponding 'Font' parameters, except that +% they apply to line widths instead of font sizes. +% +% Style Map Parameter: +% 'LineStyleMap' one of [], 'bw', or a function name or handle +% specifies how to map line colors to styles. An empty +% style map means styles are not changed. The style map +% 'bw' is a built-in mapping that maps lines with the same +% color to the same style and otherwise cycles through the +% available styles. A user-specified map is a function +% that takes as input a cell array of line objects and +% outputs a cell array of line style strings. The default +% map is []. +% +% Examples: +% exportfig(gcf,'fig1.eps','height',3); +% Exports the current figure to the file named 'fig1.eps' with +% a height of 3 inches (assuming the figure's PaperUnits is +% inches) and an aspect ratio the same as the figure's aspect +% ratio on screen. +% +% opts = struct('FontMode','fixed','FontSize',10,'height',3); +% exportfig(gcf, 'fig2.eps', opts, 'height', 5); +% Exports the current figure to 'fig2.eps' with all +% text in 10 point fonts and with height 5 inches. +% +% See also PREVIEWFIG, APPLYTOFIG, RESTOREFIG, PRINT. + +% Copyright 2000 Ben Hinkle +% Email bug reports and comments to bhinkle@mathworks.com + +if (nargin < 2) + error('Too few input arguments'); +end + +% exportfig(H, filename, [options,] ...) +H = varargin{1}; +if ~LocalIsHG(H,'figure') + error('First argument must be a handle to a figure.'); +end +filename = varargin{2}; +if ~ischar(filename) + error('Second argument must be a string.'); +end +paramPairs = {varargin{3:end}}; +if nargin > 2 + if isstruct(paramPairs{1}) + pcell = LocalToCell(paramPairs{1}); + paramPairs = {pcell{:}, paramPairs{2:end}}; + end +end +verstr = version; +majorver = str2num(verstr(1)); +defaults = []; +if majorver > 5 + if ispref('exportfig','defaults') + defaults = getpref('exportfig','defaults'); + end +elseif exist('getappdata') + defaults = getappdata(0,'exportfigdefaults'); +end +if ~isempty(defaults) + dcell = LocalToCell(defaults); + paramPairs = {dcell{:}, paramPairs{:}}; +end + +% Do some validity checking on param-value pairs +if (rem(length(paramPairs),2) ~= 0) + error(['Invalid input syntax. Optional parameters and values' ... + ' must be in pairs.']); +end + +auto.format = 'eps'; +auto.preview = 'none'; +auto.width = -1; +auto.height = -1; +auto.color = 'rgb'; +auto.defaultfontsize=10; +auto.fontsize = -1; +auto.fontmode='scaled'; +auto.fontmin = 8; +auto.fontmax = 60; +auto.defaultlinewidth = 1.0; +auto.linewidth = -1; +auto.linemode=[]; +auto.linemin = 0.5; +auto.linemax = 100; +auto.fontencoding = 'latin1'; +auto.renderer = []; +auto.resolution = []; +auto.stylemap = []; +auto.applystyle = 0; +auto.refobj = -1; +auto.bounds = 'tight'; +explicitbounds = 0; +auto.lockaxes = 1; +auto.separatetext = 0; +opts = auto; + +% Process param-value pairs +args = {}; +for k = 1:2:length(paramPairs) + param = lower(paramPairs{k}); + if ~ischar(param) + error('Optional parameter names must be strings'); + end + value = paramPairs{k+1}; + + switch (param) + case 'format' + opts.format = LocalCheckAuto(lower(value),auto.format); + if strcmp(opts.format,'preview') + error(['Format ''preview'' no longer supported. Use PREVIEWFIG' ... + ' instead.']); + end + case 'preview' + opts.preview = LocalCheckAuto(lower(value),auto.preview); + if ~strcmp(opts.preview,{'none','tiff'}) + error('Preview must be ''none'' or ''tiff''.'); + end + case 'width' + opts.width = LocalToNum(value, auto.width); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.width) + error('Width must be a numeric scalar > 0'); + end + end + case 'height' + opts.height = LocalToNum(value, auto.height); + if ~ischar(value) | ~strcmp(value,'auto') + if(~LocalIsPositiveScalar(opts.height)) + error('Height must be a numeric scalar > 0'); + end + end + case 'color' + opts.color = LocalCheckAuto(lower(value),auto.color); + if ~strcmp(opts.color,{'bw','gray','rgb','cmyk'}) + error('Color must be ''bw'', ''gray'',''rgb'' or ''cmyk''.'); + end + case 'fontmode' + opts.fontmode = LocalCheckAuto(lower(value),auto.fontmode); + if ~strcmp(opts.fontmode,{'scaled','fixed'}) + error('FontMode must be ''scaled'' or ''fixed''.'); + end + case 'fontsize' + opts.fontsize = LocalToNum(value,auto.fontsize); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.fontsize) + error('FontSize must be a numeric scalar > 0'); + end + end + case 'defaultfixedfontsize' + opts.defaultfontsize = LocalToNum(value,auto.defaultfontsize); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.defaultfontsize) + error('DefaultFixedFontSize must be a numeric scalar > 0'); + end + end + case 'fontsizemin' + opts.fontmin = LocalToNum(value,auto.fontmin); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.fontmin) + error('FontSizeMin must be a numeric scalar > 0'); + end + end + case 'fontsizemax' + opts.fontmax = LocalToNum(value,auto.fontmax); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.fontmax) + error('FontSizeMax must be a numeric scalar > 0'); + end + end + case 'fontencoding' + opts.fontencoding = LocalCheckAuto(lower(value),auto.fontencoding); + if ~strcmp(opts.fontencoding,{'latin1','adobe'}) + error('FontEncoding must be ''latin1'' or ''adobe''.'); + end + case 'linemode' + opts.linemode = LocalCheckAuto(lower(value),auto.linemode); + if ~strcmp(opts.linemode,{'scaled','fixed'}) + error('LineMode must be ''scaled'' or ''fixed''.'); + end + case 'linewidth' + opts.linewidth = LocalToNum(value,auto.linewidth); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.linewidth) + error('LineWidth must be a numeric scalar > 0'); + end + end + case 'defaultfixedlinewidth' + opts.defaultlinewidth = LocalToNum(value,auto.defaultlinewidth); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.defaultlinewidth) + error(['DefaultFixedLineWidth must be a numeric scalar >' ... + ' 0']); + end + end + case 'linewidthmin' + opts.linemin = LocalToNum(value,auto.linemin); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.linemin) + error('LineWidthMin must be a numeric scalar > 0'); + end + end + case 'linewidthmax' + opts.linemax = LocalToNum(value,auto.linemax); + if ~ischar(value) | ~strcmp(value,'auto') + if ~LocalIsPositiveScalar(opts.linemax) + error('LineWidthMax must be a numeric scalar > 0'); + end + end + case 'linestylemap' + opts.stylemap = LocalCheckAuto(value,auto.stylemap); + case 'renderer' + opts.renderer = LocalCheckAuto(lower(value),auto.renderer); + if ~ischar(value) | ~strcmp(value,'auto') + if ~strcmp(opts.renderer,{'painters','zbuffer','opengl'}) + error(['Renderer must be ''painters'', ''zbuffer'' or' ... + ' ''opengl''.']); + end + end + case 'resolution' + opts.resolution = LocalToNum(value,auto.resolution); + if ~ischar(value) | ~strcmp(value,'auto') + if ~(isnumeric(value) & (prod(size(value)) == 1) & (value >= 0)); + error('Resolution must be a numeric scalar >= 0'); + end + end + case 'applystyle' % means to apply the options and not export + opts.applystyle = 1; + case 'reference' + if ischar(value) + if strcmp(value,'auto') + opts.refobj = auto.refobj; + else + opts.refobj = eval(value); + end + else + opts.refobj = value; + end + if ~LocalIsHG(opts.refobj,'axes') + error('Reference object must evaluate to an axes handle.'); + end + case 'bounds' + opts.bounds = LocalCheckAuto(lower(value),auto.bounds); + explicitbounds = 1; + if ~strcmp(opts.bounds,{'tight','loose'}) + error('Bounds must be ''tight'' or ''loose''.'); + end + case 'lockaxes' + opts.lockaxes = LocalToNum(value,auto.lockaxes); + case 'separatetext' + opts.separatetext = LocalToNum(value,auto.separatetext); + otherwise + error(['Unrecognized option ' param '.']); + end +end + +% make sure figure is up-to-date +drawnow; + +allLines = findall(H, 'type', 'line'); +allText = findall(H, 'type', 'text'); +allAxes = findall(H, 'type', 'axes'); +allImages = findall(H, 'type', 'image'); +allLights = findall(H, 'type', 'light'); +allPatch = findall(H, 'type', 'patch'); +allSurf = findall(H, 'type', 'surface'); +allRect = findall(H, 'type', 'rectangle'); +allFont = [allText; allAxes]; +allColor = [allLines; allText; allAxes; allLights]; +allMarker = [allLines; allPatch; allSurf]; +allEdge = [allPatch; allSurf]; +allCData = [allImages; allPatch; allSurf]; + +old.objs = {}; +old.prop = {}; +old.values = {}; + +% Process format +if strncmp(opts.format,'eps',3) & ~strcmp(opts.preview,'none') + args = {args{:}, ['-' opts.preview]}; +end + +hadError = 0; +oldwarn = warning; +try + + % lock axes limits, ticks and labels if requested + if opts.lockaxes + old = LocalManualAxesMode(old, allAxes, 'TickMode'); + old = LocalManualAxesMode(old, allAxes, 'TickLabelMode'); + old = LocalManualAxesMode(old, allAxes, 'LimMode'); + end + + % Process size parameters + figurePaperUnits = get(H, 'PaperUnits'); + oldFigureUnits = get(H, 'Units'); + oldFigPos = get(H,'Position'); + set(H, 'Units', figurePaperUnits); + figPos = get(H,'Position'); + refsize = figPos(3:4); + if opts.refobj ~= -1 + oldUnits = get(opts.refobj, 'Units'); + set(opts.refobj, 'Units', figurePaperUnits); + r = get(opts.refobj, 'Position'); + refsize = r(3:4); + set(opts.refobj, 'Units', oldUnits); + end + aspectRatio = refsize(1)/refsize(2); + if (opts.width == -1) & (opts.height == -1) + opts.width = refsize(1); + opts.height = refsize(2); + elseif (opts.width == -1) + opts.width = opts.height * aspectRatio; + elseif (opts.height == -1) + opts.height = opts.width / aspectRatio; + end + wscale = opts.width/refsize(1); + hscale = opts.height/refsize(2); + sizescale = min(wscale,hscale); + old = LocalPushOldData(old,H,'PaperPositionMode', ... + get(H,'PaperPositionMode')); + set(H, 'PaperPositionMode', 'auto'); + newPos = [figPos(1) figPos(2)+figPos(4)*(1-hscale) ... + wscale*figPos(3) hscale*figPos(4)]; + set(H, 'Position', newPos); + set(H, 'Units', oldFigureUnits); + + % process line-style map + if ~isempty(opts.stylemap) & ~isempty(allLines) + oldlstyle = LocalGetAsCell(allLines,'LineStyle'); + old = LocalPushOldData(old, allLines, {'LineStyle'}, ... + oldlstyle); + newlstyle = oldlstyle; + if ischar(opts.stylemap) & strcmpi(opts.stylemap,'bw') + newlstyle = LocalMapColorToStyle(allLines); + else + try + newlstyle = feval(opts.stylemap,allLines); + catch + warning(['Skipping stylemap. ' lasterr]); + end + end + set(allLines,{'LineStyle'},newlstyle); + end + + % Process rendering parameters + switch (opts.color) + case {'bw', 'gray'} + if ~strcmp(opts.color,'bw') & strncmp(opts.format,'eps',3) + opts.format = [opts.format 'c']; + end + args = {args{:}, ['-d' opts.format]}; + + %compute and set gray colormap + oldcmap = get(H,'Colormap'); + newgrays = 0.30*oldcmap(:,1) + 0.59*oldcmap(:,2) + 0.11*oldcmap(:,3); + newcmap = [newgrays newgrays newgrays]; + old = LocalPushOldData(old, H, 'Colormap', oldcmap); + set(H, 'Colormap', newcmap); + + %compute and set ColorSpec and CData properties + old = LocalUpdateColors(allColor, 'color', old); + old = LocalUpdateColors(allAxes, 'xcolor', old); + old = LocalUpdateColors(allAxes, 'ycolor', old); + old = LocalUpdateColors(allAxes, 'zcolor', old); + old = LocalUpdateColors(allMarker, 'MarkerEdgeColor', old); + old = LocalUpdateColors(allMarker, 'MarkerFaceColor', old); + old = LocalUpdateColors(allEdge, 'EdgeColor', old); + old = LocalUpdateColors(allEdge, 'FaceColor', old); + old = LocalUpdateColors(allCData, 'CData', old); + + case {'rgb','cmyk'} + if strncmp(opts.format,'eps',3) + opts.format = [opts.format 'c']; + args = {args{:}, ['-d' opts.format]}; + if strcmp(opts.color,'cmyk') + args = {args{:}, '-cmyk'}; + end + else + args = {args{:}, ['-d' opts.format]}; + end + otherwise + error('Invalid Color parameter'); + end + if (~isempty(opts.renderer)) + args = {args{:}, ['-' opts.renderer]}; + end + if (~isempty(opts.resolution)) | ~strncmp(opts.format,'eps',3) + if isempty(opts.resolution) + opts.resolution = 0; + end + args = {args{:}, ['-r' int2str(opts.resolution)]}; + end + + % Process font parameters + if ~isempty(opts.fontmode) + oldfonts = LocalGetAsCell(allFont,'FontSize'); + oldfontunits = LocalGetAsCell(allFont,'FontUnits'); + set(allFont,'FontUnits','points'); + switch (opts.fontmode) + case 'fixed' + if (opts.fontsize == -1) + set(allFont,'FontSize',opts.defaultfontsize); + else + set(allFont,'FontSize',opts.fontsize); + end + case 'scaled' + if (opts.fontsize == -1) + scale = sizescale; + else + scale = opts.fontsize; + end + newfonts = LocalScale(oldfonts,scale,opts.fontmin,opts.fontmax); + set(allFont,{'FontSize'},newfonts); + otherwise + error('Invalid FontMode parameter'); + end + old = LocalPushOldData(old, allFont, {'FontSize'}, oldfonts); + old = LocalPushOldData(old, allFont, {'FontUnits'}, oldfontunits); + end + if strcmp(opts.fontencoding,'adobe') & strncmp(opts.format,'eps',3) + args = {args{:}, '-adobecset'}; + end + + % Process line parameters + if ~isempty(opts.linemode) + oldlines = LocalGetAsCell(allMarker,'LineWidth'); + old = LocalPushOldData(old, allMarker, {'LineWidth'}, oldlines); + switch (opts.linemode) + case 'fixed' + if (opts.linewidth == -1) + set(allMarker,'LineWidth',opts.defaultlinewidth); + else + set(allMarker,'LineWidth',opts.linewidth); + end + case 'scaled' + if (opts.linewidth == -1) + scale = sizescale; + else + scale = opts.linewidth; + end + newlines = LocalScale(oldlines, scale, opts.linemin, opts.linemax); + set(allMarker,{'LineWidth'},newlines); + end + end + + % adjust figure bounds to surround axes + if strcmp(opts.bounds,'tight') + if (~strncmp(opts.format,'eps',3) & LocalHas3DPlot(allAxes)) | ... + (strncmp(opts.format,'eps',3) & opts.separatetext) + if (explicitbounds == 1) + warning(['Cannot compute ''tight'' bounds. Using ''loose''' ... + ' bounds.']); + end + opts.bounds = 'loose'; + end + end + warning('off'); + if ~isempty(allAxes) + if strncmp(opts.format,'eps',3) + if strcmp(opts.bounds,'loose') + args = {args{:}, '-loose'}; + end + old = LocalPushOldData(old,H,'Position', oldFigPos); + elseif strcmp(opts.bounds,'tight') + oldaunits = LocalGetAsCell(allAxes,'Units'); + oldapos = LocalGetAsCell(allAxes,'Position'); + oldtunits = LocalGetAsCell(allText,'units'); + oldtpos = LocalGetAsCell(allText,'Position'); + set(allAxes,'units','points'); + apos = LocalGetAsCell(allAxes,'Position'); + oldunits = get(H,'Units'); + set(H,'units','points'); + origfr = get(H,'position'); + fr = []; + for k=1:length(allAxes) + if ~strcmpi(get(allAxes(k),'Tag'),'legend') + axesR = apos{k}; + r = LocalAxesTightBoundingBox(axesR, allAxes(k)); + r(1:2) = r(1:2) + axesR(1:2); + fr = LocalUnionRect(fr,r); + end + end + if isempty(fr) + fr = [0 0 origfr(3:4)]; + end + for k=1:length(allAxes) + ax = allAxes(k); + r = apos{k}; + r(1:2) = r(1:2) - fr(1:2); + set(ax,'Position',r); + end + old = LocalPushOldData(old, allAxes, {'Position'}, oldapos); + old = LocalPushOldData(old, allText, {'Position'}, oldtpos); + old = LocalPushOldData(old, allText, {'Units'}, oldtunits); + old = LocalPushOldData(old, allAxes, {'Units'}, oldaunits); + old = LocalPushOldData(old, H, 'Position', oldFigPos); + old = LocalPushOldData(old, H, 'Units', oldFigureUnits); + r = [origfr(1) origfr(2)+origfr(4)-fr(4) fr(3:4)]; + set(H,'Position',r); + else + args = {args{:}, '-loose'}; + old = LocalPushOldData(old,H,'Position', oldFigPos); + end + end + + % Process text in a separate file if needed + if opts.separatetext & ~opts.applystyle + % First hide all text and export + oldtvis = LocalGetAsCell(allText,'visible'); + set(allText,'visible','off'); + oldax = LocalGetAsCell(allAxes,'XTickLabel',1); + olday = LocalGetAsCell(allAxes,'YTickLabel',1); + oldaz = LocalGetAsCell(allAxes,'ZTickLabel',1); + null = cell(length(oldax),1); + [null{:}] = deal([]); + set(allAxes,{'XTickLabel'},null); + set(allAxes,{'YTickLabel'},null); + set(allAxes,{'ZTickLabel'},null); + print(H, filename, args{:}); + set(allText,{'Visible'},oldtvis); + set(allAxes,{'XTickLabel'},oldax); + set(allAxes,{'YTickLabel'},olday); + set(allAxes,{'ZTickLabel'},oldaz); + % Now hide all non-text and export as eps in painters + [path, name, ext] = fileparts(filename); + tfile = fullfile(path,[name '_t.eps']); + tfile2 = fullfile(path,[name '_t2.eps']); + foundRenderer = 0; + for k=1:length(args) + if strncmp('-d',args{k},2) + args{k} = '-deps'; + elseif strncmp('-zbuffer',args{k},8) | ... + strncmp('-opengl', args{k},6) + args{k} = '-painters'; + foundRenderer = 1; + end + end + if ~foundRenderer + args = {args{:}, '-painters'}; + end + allNonText = [allLines; allLights; allPatch; ... + allImages; allSurf; allRect]; + oldvis = LocalGetAsCell(allNonText,'visible'); + oldc = LocalGetAsCell(allAxes,'color'); + oldaxg = LocalGetAsCell(allAxes,'XGrid'); + oldayg = LocalGetAsCell(allAxes,'YGrid'); + oldazg = LocalGetAsCell(allAxes,'ZGrid'); + [null{:}] = deal('off'); + set(allAxes,{'XGrid'},null); + set(allAxes,{'YGrid'},null); + set(allAxes,{'ZGrid'},null); + set(allNonText,'Visible','off'); + set(allAxes,'Color','none'); + print(H, tfile2, args{:}); + set(allNonText,{'Visible'},oldvis); + set(allAxes,{'Color'},oldc); + set(allAxes,{'XGrid'},oldaxg); + set(allAxes,{'YGrid'},oldayg); + set(allAxes,{'ZGrid'},oldazg); + %hack up the postscript file + fid1 = fopen(tfile,'w'); + fid2 = fopen(tfile2,'r'); + line = fgetl(fid2); + while ischar(line) + if strncmp(line,'%%Title',7) + fprintf(fid1,'%s\n',['%%Title: ', tfile]); + elseif (length(line) < 3) + fprintf(fid1,'%s\n',line); + elseif ~strcmp(line(end-2:end),' PR') & ... + ~strcmp(line(end-1:end),' L') + fprintf(fid1,'%s\n',line); + end + line = fgetl(fid2); + end + fclose(fid1); + fclose(fid2); + delete(tfile2); + + elseif ~opts.applystyle + drawnow; + print(H, filename, args{:}); + end + warning(oldwarn); + +catch + warning(oldwarn); + hadError = 1; +end + +% Restore figure settings +if opts.applystyle + varargout{1} = old; +else + for n=1:length(old.objs) + if ~iscell(old.values{n}) & iscell(old.prop{n}) + old.values{n} = {old.values{n}}; + end + set(old.objs{n}, old.prop{n}, old.values{n}); + end +end + +if hadError + error(deblank(lasterr)); +end + +% +% Local Functions +% + +function outData = LocalPushOldData(inData, objs, prop, values) +outData.objs = {objs, inData.objs{:}}; +outData.prop = {prop, inData.prop{:}}; +outData.values = {values, inData.values{:}}; + +function cellArray = LocalGetAsCell(fig,prop,allowemptycell); +cellArray = get(fig,prop); +if nargin < 3 + allowemptycell = 0; +end +if ~iscell(cellArray) & (allowemptycell | ~isempty(cellArray)) + cellArray = {cellArray}; +end + +function newArray = LocalScale(inArray, scale, minv, maxv) +n = length(inArray); +newArray = cell(n,1); +for k=1:n + newArray{k} = min(maxv,max(minv,scale*inArray{k}(1))); +end + +function gray = LocalMapToGray1(color) +gray = color; +if ischar(color) + switch color(1) + case 'y' + color = [1 1 0]; + case 'm' + color = [1 0 1]; + case 'c' + color = [0 1 1]; + case 'r' + color = [1 0 0]; + case 'g' + color = [0 1 0]; + case 'b' + color = [0 0 1]; + case 'w' + color = [1 1 1]; + case 'k' + color = [0 0 0]; + end +end +if ~ischar(color) + gray = 0.30*color(1) + 0.59*color(2) + 0.11*color(3); +end + +function newArray = LocalMapToGray(inArray); +n = length(inArray); +newArray = cell(n,1); +for k=1:n + color = inArray{k}; + if ~isempty(color) + color = LocalMapToGray1(color); + end + if isempty(color) | ischar(color) + newArray{k} = color; + else + newArray{k} = [color color color]; + end +end + +function newArray = LocalMapColorToStyle(inArray); +inArray = LocalGetAsCell(inArray,'Color'); +n = length(inArray); +newArray = cell(n,1); +styles = {'-','--',':','-.'}; +uniques = []; +nstyles = length(styles); +for k=1:n + gray = LocalMapToGray1(inArray{k}); + if isempty(gray) | ischar(gray) | gray < .05 + newArray{k} = '-'; + else + if ~isempty(uniques) & any(gray == uniques) + ind = find(gray==uniques); + else + uniques = [uniques gray]; + ind = length(uniques); + end + newArray{k} = styles{mod(ind-1,nstyles)+1}; + end +end + +function newArray = LocalMapCData(inArray); +n = length(inArray); +newArray = cell(n,1); +for k=1:n + color = inArray{k}; + if (ndims(color) == 3) & isa(color,'double') + gray = 0.30*color(:,:,1) + 0.59*color(:,:,2) + 0.11*color(:,:,3); + color(:,:,1) = gray; + color(:,:,2) = gray; + color(:,:,3) = gray; + end + newArray{k} = color; +end + +function outData = LocalUpdateColors(inArray, prop, inData) +value = LocalGetAsCell(inArray,prop); +outData.objs = {inData.objs{:}, inArray}; +outData.prop = {inData.prop{:}, {prop}}; +outData.values = {inData.values{:}, value}; +if (~isempty(value)) + if strcmp(prop,'CData') + value = LocalMapCData(value); + else + value = LocalMapToGray(value); + end + set(inArray,{prop},value); +end + +function bool = LocalIsPositiveScalar(value) +bool = isnumeric(value) & ... + prod(size(value)) == 1 & ... + value > 0; + +function value = LocalToNum(value,auto) +if ischar(value) + if strcmp(value,'auto') + value = auto; + else + value = str2num(value); + end +end + +%convert a struct to {field1,val1,field2,val2,...} +function c = LocalToCell(s) +f = fieldnames(s); +v = struct2cell(s); +opts = cell(2,length(f)); +opts(1,:) = f; +opts(2,:) = v; +c = {opts{:}}; + +function c = LocalIsHG(obj,hgtype) +c = 0; +if (length(obj) == 1) & ishandle(obj) + c = strcmp(get(obj,'type'),hgtype); +end + +function c = LocalHas3DPlot(a) +zticks = LocalGetAsCell(a,'ZTickLabel'); +c = 0; +for k=1:length(zticks) + if ~isempty(zticks{k}) + c = 1; + return; + end +end + +function r = LocalUnionRect(r1,r2) +if isempty(r1) + r = r2; +elseif isempty(r2) + r = r1; +elseif max(r2(3:4)) > 0 + left = min(r1(1),r2(1)); + bot = min(r1(2),r2(2)); + right = max(r1(1)+r1(3),r2(1)+r2(3)); + top = max(r1(2)+r1(4),r2(2)+r2(4)); + r = [left bot right-left top-bot]; +else + r = r1; +end + +function c = LocalLabelsMatchTicks(labs,ticks) +c = 0; +try + t1 = num2str(ticks(1)); + n = length(ticks); + tend = num2str(ticks(n)); + c = strncmp(labs(1),t1,length(labs(1))) & ... + strncmp(labs(n),tend,length(labs(n))); +end + +function r = LocalAxesTightBoundingBox(axesR, a) +r = []; +atext = findall(a,'type','text','visible','on'); +if ~isempty(atext) + set(atext,'units','points'); + res=LocalGetAsCell(atext,'extent'); + for n=1:length(atext) + r = LocalUnionRect(r,res{n}); + end +end +if strcmp(get(a,'visible'),'on') + r = LocalUnionRect(r,[0 0 axesR(3:4)]); + oldunits = get(a,'fontunits'); + set(a,'fontunits','points'); + label = text(0,0,'','parent',a,... + 'units','points',... + 'fontsize',get(a,'fontsize'),... + 'fontname',get(a,'fontname'),... + 'fontweight',get(a,'fontweight'),... + 'fontangle',get(a,'fontangle'),... + 'visible','off'); + fs = get(a,'fontsize'); + + % handle y axis tick labels + ry = [0 -fs/2 0 axesR(4)+fs]; + ylabs = get(a,'yticklabels'); + yticks = get(a,'ytick'); + maxw = 0; + if ~isempty(ylabs) + for n=1:size(ylabs,1) + set(label,'string',ylabs(n,:)); + ext = get(label,'extent'); + maxw = max(maxw,ext(3)); + end + if ~LocalLabelsMatchTicks(ylabs,yticks) & ... + strcmp(get(a,'xaxislocation'),'bottom') + ry(4) = ry(4) + 1.5*ext(4); + end + if strcmp(get(a,'yaxislocation'),'left') + ry(1) = -(maxw+5); + else + ry(1) = axesR(3); + end + ry(3) = maxw+5; + r = LocalUnionRect(r,ry); + end + + % handle x axis tick labels + rx = [0 0 0 fs+5]; + xlabs = get(a,'xticklabels'); + xticks = get(a,'xtick'); + if ~isempty(xlabs) + if strcmp(get(a,'xaxislocation'),'bottom') + rx(2) = -(fs+5); + if ~LocalLabelsMatchTicks(xlabs,xticks); + rx(4) = rx(4) + 2*fs; + rx(2) = rx(2) - 2*fs; + end + else + rx(2) = axesR(4); + % exponent is still below axes + if ~LocalLabelsMatchTicks(xlabs,xticks); + rx(4) = rx(4) + axesR(4) + 2*fs; + rx(2) = -2*fs; + end + end + set(label,'string',xlabs(1,:)); + ext1 = get(label,'extent'); + rx(1) = -ext1(3)/2; + set(label,'string',xlabs(size(xlabs,1),:)); + ext2 = get(label,'extent'); + rx(3) = axesR(3) + (ext2(3) + ext1(3))/2; + r = LocalUnionRect(r,rx); + end + set(a,'fontunits',oldunits); + delete(label); +end + +function c = LocalManualAxesMode(old, allAxes, base) +xs = ['X' base]; +ys = ['Y' base]; +zs = ['Z' base]; +oldXMode = LocalGetAsCell(allAxes,xs); +oldYMode = LocalGetAsCell(allAxes,ys); +oldZMode = LocalGetAsCell(allAxes,zs); +old = LocalPushOldData(old, allAxes, {xs}, oldXMode); +old = LocalPushOldData(old, allAxes, {ys}, oldYMode); +old = LocalPushOldData(old, allAxes, {zs}, oldZMode); +set(allAxes,xs,'manual'); +set(allAxes,ys,'manual'); +set(allAxes,zs,'manual'); +c = old; + +function val = LocalCheckAuto(val, auto) +if ischar(val) & strcmp(val,'auto') + val = auto; +end diff --git a/sourcecodes/bnt-master/SLP/misc/find_nodes_in_undirected_component.m b/sourcecodes/bnt-master/SLP/misc/find_nodes_in_undirected_component.m new file mode 100644 index 00000000..7ba14eca --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/find_nodes_in_undirected_component.m @@ -0,0 +1,18 @@ +function unprocessed = find_unprocessed(dag) +% unprocessed = find_unprocessed(dag) +% +% francois.olivier.c.h@gmail.com + +undirected_edges = dag.*dag'; +[unprocessed, tmp] = find(undirected_edges) +unprocessed = unique(unprocessed); +% N = size(dag,1); +% unprocessed = []; +% for i=1:(N-1) +% for j=(i+1):N +% if dag(i,j)==1 & dag(j,i)==1 +% unprocessed = [unprocessed,i,j]; +% end +% end +% end +% unprocessed = unique(unprocessed); diff --git a/sourcecodes/bnt-master/SLP/misc/foptions.m b/sourcecodes/bnt-master/SLP/misc/foptions.m new file mode 100644 index 00000000..f2976bb9 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/foptions.m @@ -0,0 +1,11 @@ +function opt_vect = foptions() +% FOPTIONS Sets default parameters for optimisation routines +% For compatibility with MATLAB's foptions() +% +% Copyright (c) Dharmesh Maniyar, Ian T. Nabney (2004) + +opt_vect = zeros(1, 18); +opt_vect(2:3) = 1e-4; +opt_vect(4) = 1e-6; +opt_vect(16) = 1e-8; +opt_vect(17) = 0.1; \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/gener_MAR_net.m b/sourcecodes/bnt-master/SLP/misc/gener_MAR_net.m new file mode 100644 index 00000000..64740076 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/gener_MAR_net.m @@ -0,0 +1,104 @@ +function bnet_miss = gener_MAR_net(bnet_orig, base_proba) +% function bnet_miss = gener_MAR_net(bnet_orig, base_proba) +% +% bnet_orig : a bnet +% base_proba : a probability for value to be missing +% +% bnet_miss : a bnet that could be used in gener_data_from_bnet_miss function +% to generate incomplete MAR dataset +% +% Francois.Olivier.C.H@gmail.com + +%%%%%%%%%%%% INIT +if nargin<2, error('Not enougth arguments'); end + +% création du réseau +dag = bnet_orig.dag; +N = size(dag,2); +ns = bnet_orig.node_sizes; + + ns_miss = zeros(1,3*N); + ns_miss(1:N) = ns; + ns_miss(N+1:2*N) = 2*ones(1,N); % 1= node i-N present, 2= node i-N missing + ns_miss(2*N+1:3*N) = ns+1; % 1:ns, absent + + dag_miss = zeros(3*N,3*N); + dag_miss(1:N,1:N) = dag; +% dag_miss(2*N+1:3*N,N+1:2*N)=mk_rnd_dag(N,N-ceil(rand*N/2)); + lim = 1+(rand>.4)+(rand>.65)+(rand>.9); + dag_miss(2*N+1:3*N,N+1:2*N)=mk_rnd_dag(N,lim); + + for i=1:N, dag_miss(i,2*N+i)=1; dag_miss(N+i,2*N+i)=1; dag_miss(2*N+i,i)=0; dag_miss(2*N+i,N+i)=0; end + + bnet_miss = mk_bnet(dag_miss, ns_miss); + CPT = CPT_from_bnet(bnet_orig, 0); + for i=1:N + bnet_miss.CPD{i} = tabular_CPD (bnet_miss, i, CPT{i}); + end + + % CPD of nodes M + for i=1:N, + if find(bnet_miss.order==i)<find(bnet_miss.order==N+i), + CPT_M=[]; + for j=1:ns_miss(i), for l=1:ns_miss(N+i), for k=1:ns_miss(2*N+i), + CPT_M=[CPT_M (((j==k)&(l==1))|((k==ns_miss(2*N+i))&(l==2)))]; + end, end, end + else + CPT_M=[]; + for k=1:ns_miss(2*N+i), for l=1:ns_miss(N+i), for j=1:ns_miss(i), + CPT_M=[CPT_M (((j==k)&(l==1))|((k==ns_miss(2*N+i))&(l==2)))]; + end, end, end + end + bnet_miss.CPD{2*N+i} = tabular_CPD (bnet_miss, 2*N+i, CPT_M); + end + +%%%%%%%%%%%% Base probability of missing value +p = base_proba; + for i=1:N + fam = find(dag_miss(:,N+i)==1)'; + semisize = prod(ns_miss(fam)); % as node N+i is binary to say i is present or missing + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-p; + CPT(semisize+1:2*semisize) = p; + bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); + end + + BETA = gener_discrete_dist(N, base_proba); + + order=[]; + missdagtmp = bnet_miss.dag; %(N+1:2*N,N+1:2*N); + unprocessed = 1:N; + while ~isempty(unprocessed) + npar=[]; + for i=N+1:2*N, npar(end+1)=length(parents(missdagtmp,i));end, % to be verifie from here + [npar, ord] = sort(npar); + while ~ismember(ord(1),unprocessed) + ord=ord(2:end); + end + order = [order, ord(1)]; + missdagtmp(ord(1),:)=0; + unprocessed = mysetdiff(unprocessed,ord(1)); + end + +%%%%%%%%%%%% Update CPT with MCAR process +for i=1:length(BETA) + fam_miss = find(dag_miss(:,N+order(i))==1)'; + p=BETA(i); + semisize = prod(ns_miss(fam_miss)); + + if isempty(fam_miss), + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-p; + CPT(semisize+1:2*semisize) = p; + bnet_miss.CPD{N+order(i)} = tabular_CPD (bnet_miss, N+order(i), CPT); + else + + MUi1k = gener_discrete_dist(semisize, p); + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-MUi1k; + CPT(semisize+1:2*semisize) = MUi1k; + bnet_miss.CPD{N+order(i)} = tabular_CPD (bnet_miss, N+order(i), CPT); + + end + +end diff --git a/sourcecodes/bnt-master/SLP/misc/gener_MCAR_net.m b/sourcecodes/bnt-master/SLP/misc/gener_MCAR_net.m new file mode 100644 index 00000000..5e852a23 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/gener_MCAR_net.m @@ -0,0 +1,165 @@ +function bnet_miss = gener_MCAR_net(bnet_orig, base_proba) +% function bnet_miss = gener_MCAR_net(bnet_orig, base_proba) +% +% bnet_orig : a bnet +% base_proba : a probability for value to be missing +% +% bnet_miss : a bnet that could be used in gener_data_from_bnet_miss function +% to generate incomplete MCAR dataset +% +% Francois.Olivier.C.H@gmail.com + +%%%%%%%%%%%% INIT +%bnet_miss = gener_MCAR_net(bnet_orig, base_proba, upd, bnet_biss) +%if nargin<4, upd = 0; else upd = 1; end +%if nargin<3, manual = 0; end +if nargin<2, error('Not enougth arguments'); end + +% création du réseau +dag = bnet_orig.dag; +N = size(dag,2); +ns = bnet_orig.node_sizes; + +%if nargin<4, + ns_miss = zeros(1,3*N); + ns_miss(1:N) = ns; + ns_miss(N+1:2*N) = 2*ones(1,N); % 1= node i-N present, 2= node i-N missing + ns_miss(2*N+1:3*N) = ns+1; % 1:ns, absent + + dag_miss = zeros(3*N,3*N); + dag_miss(1:N,1:N) = dag; + dag_miss(N+1:2*N,N+1:2*N)=mk_rnd_dag(N,N-ceil(rand*N/2)); %dag_miss(N+1:2*N,N+1:2*N)=mk_rnd_dag(N,2); + for i=1:N, dag_miss(i,2*N+i)=1; dag_miss(N+i,2*N+i)=1; end + + bnet_miss = mk_bnet(dag_miss, ns_miss); + CPT = CPT_from_bnet(bnet_orig, 0); + for i=1:N + bnet_miss.CPD{i} = tabular_CPD (bnet_miss, i, CPT{i}); + end + + % CPD of nodes M + for i=1:N, + if find(bnet_miss.order==i)<find(bnet_miss.order==N+i), + CPT_M=[]; + for j=1:ns_miss(i), for l=1:ns_miss(N+i), for k=1:ns_miss(2*N+i), + CPT_M=[CPT_M (((j==k)&(l==1))|((k==ns_miss(2*N+i))&(l==2)))]; + end, end, end + else + CPT_M=[]; + for k=1:ns_miss(2*N+i), for l=1:ns_miss(N+i), for j=1:ns_miss(i), + CPT_M=[CPT_M (((j==k)&(l==1))|((k==ns_miss(2*N+i))&(l==2)))]; + end, end, end + end + bnet_miss.CPD{2*N+i} = tabular_CPD (bnet_miss, 2*N+i, CPT_M); + end +% else +% ns_miss = bnet_miss.node_sizes; +% dag_miss = bnet_miss.dag; +% end + +%%%%%%%%%%%% Base probability of missing value +p = base_proba; + for i=1:N + fam = find(dag_miss(:,N+i)==1)'; + semisize = prod(ns_miss(fam)); % as node N+i is binary to say i is present or missing + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-p; + CPT(semisize+1:2*semisize) = p; + bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); + end + +% if manual, % manual generation +% +% if upd==1; +% +% fprintf('Base probability of a data to be missing is %1.4f',base_proba); +% +% b=-1; +% while ~(b==0 | b==1), b = input('Would you like to change a probability of a node to be missing (1 for yes, 0 for no) ? '); end +% +% %%%%%%%%%%%% Update CPT with MCAR process +% while b +% fprintf('Nodes are from 1 to %d. ',N); +% i=0; +% while i<1 | i>N | round(i)~=i, i = input('Which node ? '); end +% fam = find(dag(:,i)==1)'; +% fam_miss = find(dag_miss(:,N+i)==1)'; +% cas = -ones(1, length(fam_miss)+1); +% familly = [fam_miss, N+i]; +% p=-1; +% while p<0 | p>1, p = input(' - A priori probability for this node to be missing ? ');end +% semisize = prod(ns_miss(fam_miss)); +% CPT = zeros(1,2*semisize); +% CPT(1:semisize) = 1-p; +% CPT(semisize+1:2*semisize) = p; +% bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); +% b=-1; +% while ~(b==0 | b==1), b = input('Would you like to change a probability of a node to be missing (1 for yes, 0 for no) ? ');end +% end +% end +% +% +% else % automatic generation + +%%%%%%%%%%%% creating BETAs +% %% To use multiple of 5 percent in probs +% if N<=25, +% BETA = gener_problist(base_proba, N); +% else +% nboucles = floor(N/25); +% BETA = []; +% for i=1:nboucles +% BETA1 = gener_problist(base_proba, 25); +% BETA = [BETA, BETA1]; +% end +% nreste = rem(N,25); +% BETA1 = gener_problist(base_proba, nreste); +% BETA = [BETA, BETA1]; +% end + BETA = gener_discrete_dist(N, base_proba); + + order=[]; + missdagtmp = bnet_miss.dag(N+1:2*N,N+1:2*N); + unprocessed = 1:N; + while ~isempty(unprocessed) + npar=[]; + for i=1:N, npar(end+1)=length(parents(missdagtmp,i));end, + [npar, ord] = sort(npar); + while ~ismember(ord(1),unprocessed) + ord=ord(2:end); + end + order = [order, ord(1)]; + missdagtmp(ord(1),:)=0; + unprocessed = mysetdiff(unprocessed,ord(1)); + end + +%%%%%%%%%%%% Update CPT with MCAR process +for i=1:length(BETA) + fam_miss = find(dag_miss(:,N+order(i))==1)'; + p=BETA(i); + semisize = prod(ns_miss(fam_miss)); + + if isempty(fam_miss), + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-p; + CPT(semisize+1:2*semisize) = p; + bnet_miss.CPD{N+order(i)} = tabular_CPD (bnet_miss, N+order(i), CPT); + else + + %node = N+order(i) + %for k=1:semisize + % XI(k) = eval_xi(bnet_miss, N+order(i), k); + % %XI(k+semisize)=1-XI(k); + %end + %MUi1 = zeros(1,semisize); + %MUi1 = gener_mu(p, semisize, XI); + + MUi1k = gener_discrete_dist(semisize, p); + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-MUi1k; + CPT(semisize+1:2*semisize) = MUi1k; + bnet_miss.CPD{N+order(i)} = tabular_CPD (bnet_miss, N+order(i), CPT); + + end +% end +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/gener_NMAR_data.m b/sourcecodes/bnt-master/SLP/misc/gener_NMAR_data.m new file mode 100644 index 00000000..4cff44d6 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/gener_NMAR_data.m @@ -0,0 +1,290 @@ +function [data, comp_data, bnet_miss] = gener_NMAR_dataset(bnet_orig, m, bnet_miss, upd) +% [NMAR_data] = gener_NMAR_data(bnet_miss, length_of_dataset) +% +% this function takes in input a bnet that can be used to generate NMAR data. +% this bnet bnet_miss can be creating by the function gener_NMAR_bnet. +% +% - comp_data (array) is a dataset that was generate by the bnet_orig that you enter in gener_NMAR_bnet +% - NMAR_data (cell array) is the dataset compdata that was emptyied by the NMAR process encodes in bnet_miss +% +% optional : +% - bnet_miss : an old bnet_miss built by this function +% - upd==1 if you want to update the bnet_miss +% +% [data, comp_data, bnet_miss] = gener_NMAR_dataset(bnet_orig, m, bnet_miss, upd); +% +% version 0.5 : june 8th 2005, olivier.francois@insa-rouen.fr +% +% TO DO : +% - allow the combinaison of a missing state of one variable and another state of another variable to have influence +% - allow the introduction of new nodes and specify which nodes it influence and which nodes has influence on it (it will also satisfy the first task then) +% + + +% INIT +N = size(bnet_orig.dag,2); +if mod(N,2)~=0, error('The number of nodes must be even'); end +if nargin<4, upd =0; end + +% fisrt rules +if nargin<3, + l1=[]; l2=[]; lp=[]; b=-1; + while ~(b==0 | b==1), b = input('Would you like to make a node missing when another one is missing (1 for yes, 0 for no) ? '); end +else + b=-1; + if upd, while ~(b==0 | b==1), b = input('Would you like to add rules (1 for yes, 0 for no) ? '); end + l1=bnet_miss.list{1}; + l2=bnet_miss.list{2}; + lp=bnet_miss.list{3}; + else + l1=[]; l2=[]; lp=[]; b=-1; + end +end + +while b==1, + n = input('The firts node ? '); + s = input('The node that have to be missing when this one is missing ? '); + p = input('The probability of the second node to be missing ? '); + l1 = [l1, n]; l2 = [l2, s]; lp=[lp, p]; + b=-1; + while ~(b==0 | b==1), b = input('Another one (1 for yes, 0 for no) ? '); end +end +bb = length(l1); + + +if nargin>=3, + bnet_miss = gener_NMAR_bnet(upd, bnet_orig, bnet_miss); +else + bnet_miss = gener_NMAR_bnet(1, bnet_orig); + bnet_miss.list={l1; l2; lp}; +%%%%%%%%%%% SAVING FILE + ss = 1; + if nargin == 2 | upd==1, ss = input('Would you like to save the bnet of the NMAR process you have made (1 for yes) ? '); end + if ss == 1, + ddd = datestr(now); + ddd([12 15 18])='-' ; + fnout=['NMAR-bnet-' ddd '.mat']; + eval(['save ' fnout ' bnet_miss']); + fprintf(' The bnet for NMAR process was saved as : %s\n',fnout); + end +end %if nargin + +% Generation of a complete dataset +if N>9 & m>2000, disp(' ! It could take a long time...'); end +data = cell(N,m); +for l = 1:m, data(:,l) = sample_bnet(bnet_orig); end +disp('Complete data have been creating.'); + +% Generation of a NMAR dataset +miss_array = cell(2*N,m); +vide = cell(1,N); l= 1; +while l <= m, + ev(1:N) = data(:,l); ev(N+1:2*N) = vide; + miss_array(:,l) = sample_bnet(bnet_miss, 'evidence', ev); + % apply simple rule of missingness + ev2 = cell2mat(miss_array(N+1:2*N, l)); + if bb, + missl1 = myintersect(find(ev2==2), l1); + if ~isempty(missl1), + for i=1:length(l1), + if ev2(l1(i))==2, if rand<lp(i), ev2(l2(i))=2; miss_array{N+l2(i),l}=[2]; end, end + end, end, end + % verification that we have not a completly missing sample + ev2 = 3-ev2; + if prod(ev2)==1, fprintf(' - %d, one completly missing sample removed', l); else l=l+1; end + if mod(l,100)==0, fprintf('\n - %d',l); end +end +fprintf('\n'); +data = bnt_to_mat(data); comp_data = data; +miss_array = bnt_to_mat(miss_array(N+1:2*N, :)); +miss_array = 2-miss_array; +data = data.*miss_array; +data = mat_to_bnt(data, 0); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function bnet_miss = gener_NMAR_bnet(upd, bnet_orig, bnet_miss) + +%%%%%%%%%%%% INIT +dag = bnet_orig.dag; +N = size(dag,2); +ns = bnet_orig.node_sizes; +NN = 0; +ns_miss = zeros(1,2*N+NN); +ns_miss(1:N) = ns; +ns_miss(N+1:2*N) = 2*ones(1,N); % 1= node i-N present, 2= node i-N missing + dagm = zeros(2*N,2*N); + dagm(1:N,1:N) = dag; + dagm(N+1:2*N, N+1:2*N) = dag; + dagm(1:N, N+1:2*N) = dag; + for i=1:N, dagm(i, i+N)=1; end + +%%%%%%%%%%%% NEW NODES +% b=-1; +% while ~(b==0 | b==1), b = input('Would you like to add new nodes ? '); end +% +% if b==1, +% NN = -1; +% while (NN<0 | round(NN)~=NN), NN = input('How Many ? '); end +% if NN==1, fprintf('New node will be called %d\n',(2*N+1)); +% else fprintf('New node will be called %d and following numbers\n',(2*N+1)); end +% +% disp(' !!! Make sur that the dependence you will create will not create cycle in the Bnet used to create NMAR data !!! here is the current DAG'); +% dagm2 = zeros(2*N+NN, 2*N+NN); + dagm2(1:2*N, 1:2*N) = dagm; +% draw_graph(dagm2); drawnow; +% clear dagm +% +% L1={};L2={}; +% for i=1:NN +% fprintf('For the node %d, ',(2*N+i)); +% L1{i} = input('FROM which nodes will it have influence (sample [3 2]) ? '); +% L2{i} = input('ON which nodes will it have influence (sample [1 4 3]) ? '); +% ns_miss(2*N+i) = input('What is its size ? '); +% end +% +% end + +%%%%%%%%%%%% BNET CREATION +if nargin==2, + dag_miss = dagm2; + for i=1:NN, dag_miss(L1{i},2*N+i)=1; dag_miss(2*N+i, L2{i})=1; end + + bnet_miss = mk_bnet(dag_miss, ns_miss); + CPT = CPT_from_bnet(bnet_orig); + for i=1:N + bnet_miss.CPD{i} = tabular_CPD (bnet_miss, i, CPT{i}); % error with new nodes + end +elseif nargin==3, + ns_miss = bnet_miss.node_sizes; + dag_miss = bnet_miss.dag; +end +clear dagm2 +order = 1:(2*N+NN); + +%%%%%%%%%%%% Base probability of missing value +if nargin==2, b=1; else b=0; end +if b==1, +disp('Probability MUST be between 0 and 1.'); + +if nargin==2, + p=-1; + while p<0 | p>1, p = input('Base probability of a value to be missing ? '); end + for i=1:N + fam = find(dag_miss(:,N+i)==1)'; + semisize = prod(ns_miss(fam)); % as node N+i is binary to say i is present or missing + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-p; + CPT(semisize+1:2*semisize) = p; + bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); + end +end +end + +b=-1; +while ~(b==0 | b==1), b = input('Would you like to change a probability of a node to be missing (1 for yes, 0 for no) ? '); end +if b, disp(' BE CAREFULL !! New rules can overwrite old ones partialy or fully !! So the order of entries is important'); end + +%%%%%%%%%%%% Update CPT with NMAR process +while b + fprintf('Nodes are from 1 to %d. ',N); + i=0; + while i<1 | i>N | round(i)~=i, i = input('Which node ? '); end + fam = find(dag(:,i)==1)'; + fam_miss = find(dag_miss(:,N+i)==1)'; + cas = -ones(1, length(fam_miss)+1); + familly = [fam_miss, N+i]; + fprintf('States are from 1 to %d (-1 for any states, -2 to cancel). For which state of the variable %d ?', ns(i), i); + state=-3; + while state<-2 | state>ns(i) | round(state)~=state | state==0, state = input(' ');end + if isempty(fam), + if state==-1, + p=-1; + while p<0 | p>1, p = input(' - A priori probability for this node to be missing ? ');end + semisize = prod(ns_miss(fam_miss)); + CPT = zeros(1,2*semisize); + CPT(1:semisize) = 1-p; + CPT(semisize+1:2*semisize) = p; + bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); + elseif state~=-2 + cas = state; + CPT = CPT_from_bnet(bnet_miss); + CPT = CPT{N+i}; + p=-1; + while p<0 | p>1, p = input(' - A priori probability for this node to be missing in this state ? ');end + ind = subv2ind(ns_miss(familly),[cas, 1]); + CPT(ind)=1-p; + ind = subv2ind(ns_miss(familly),[cas, 2]); + CPT(ind)=p; + bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); + end + else + if state>-2, + siz=length(cas); + place = find(fam_miss==i); + cas(place) = state; + + for k = fam, + state=-3; + fprintf(' - For the parent named %d, states are from 1 to %d (-1 for any states of this parent). ',k, ns(k)); + while state<0 | state>ns(k) | round(state)~=state, state = input('Which state ? ');end + %if state==0, + % place = find(fam_miss==(fam_miss(k)+N)); + % cas(place) = 2; % Missing + %elseif state==-2, % a changer ??? + % disp(' This case is buggy, taking missing state instand to minimise influence.'); + % place = find(fam_miss==(fam_miss(k)+N)); + % cas(place) = 2; + %elseif state~=0 & state~=-2, + place = find(fam_miss==(fam_miss(k))); + cas(place) = state; % Present + if state~=-1; place = find(fam_miss==(fam_miss(k)+N)); cas(place) = 1; end + %end + end + + p=-1; + while p<0 | p>1, p = input('Probability in this case of the value to be missing ? ');end + CPT = CPT_from_bnet(bnet_miss); + CPT = CPT{N+i}; + + cas(end) = 1; % i is present + subcas_names = find(cas==-1); + if isempty(subcas_names), + ind = subv2ind(ns_miss(familly),cas); + CPT(ind) = 1-p; + else + subcas = ones(1, length(subcas_names)); + continu = 1; + while continu + cas(subcas_names) = subcas; + ind = subv2ind(ns_miss(familly),cas); + CPT(ind) = 1-p; + [subcas, continu] = next_case(subcas, ns_miss(familly(subcas_names))); + end + end + + cas(end)=2; % i is missing + if isempty(subcas_names), + ind = subv2ind(ns_miss(familly),cas); + CPT(ind) = p; + else + subcas = ones(1, length(subcas_names)); + continu = 1; + while continu + cas(subcas_names) = subcas; + ind = subv2ind(ns_miss(familly),cas); + CPT(ind) = p; + [subcas, continu] = next_case(subcas, ns_miss(familly(subcas_names))); + end + end + mass = sum(CPT, length(size(CPT))); + while length(size(mass))>2, mass = prod(mass, length(size(mass))); end + mass = prod(prod(mass)); + if mass~=1, disp('not a proba...'); end + bnet_miss.CPD{N+i} = tabular_CPD (bnet_miss, N+i, CPT); + end %if state~=-2 for the node + end %if isempty(fam), + b=-1; + while ~(b==0 | b==1), b = input('Would you like to change a probability of a node to be missing (1 for yes, 0 for no) ? ');end +end + diff --git a/sourcecodes/bnt-master/SLP/misc/gener_data_from_bnet_miss.m b/sourcecodes/bnt-master/SLP/misc/gener_data_from_bnet_miss.m new file mode 100644 index 00000000..3508daa0 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/gener_data_from_bnet_miss.m @@ -0,0 +1,103 @@ +function [data, comp_data, bnet_miss, taux, bnet_orig, notok, d] = gener_data_from_bnet_miss(bnet_miss, m, base_proba ,v, testdata) +% [data, comp_data, bnet_miss, taux, bnet_orig, notok] = gener_data_from_bnet_miss(bnet_miss, m, base_proba ,v, aretestdata) +% +% bnet_miss : see gener_[MCAR or MAR]_net function +% m : the length of the dataset +% if base_proba==0 or does not exist the x2_test will be passed +% v==1 to enter the verbose mode +% aretestdata==1 to always build the same dataset <-- rand('state',0) +% +% Francois.Olivier.C.H@gmail.com + +% Initialisation + +if nargin<5, testdata = 0; end +if nargin<4, v = 0; end + +N2 = length(bnet_miss.dag); +if mod(N2,3)~=0, error('The number of nodes must be even in bnet_miss'); end +N = length(bnet_miss.dag)/3; + +if nargin<3, base_proba=0; end +if nargin<2, error('Not enougth parameters'); end + +% CHOOSE THE TEST POWER (only affect 'notok' value) + chi2_0_1_1fd = 2.705 ; + chi2_0_05_1fd = 3.841 ; + chi2_0_01_1fd = 6.635 ; + chi2_0_001_1fd = 10.827 ; + chi2_0_0001_1fd = 15.137 ; + choice = chi2_0_001_1fd; +notok=0; +clear chi2_0_1_1fd chi2_0_05_1fd chi2_0_01_1fd chi2_0_001_1fd chi2_0_0001_1fd + +% Recovering bnet_orig + +dag = bnet_miss.dag(1:N,1:N); +bnet_orig = mk_bnet(dag, bnet_miss.node_sizes(1:N)); +CPT = CPT_from_bnet(bnet_miss, 0); +for i=1:N, bnet_orig.CPD{i} = tabular_CPD(bnet_orig, i, CPT{i}); end + +% Generation of complete data + +if testdata, rand('state',0); randn('state',0); end + +data = cell(N,m); +for l = 1:m, data(:,l) = sample_bnet(bnet_orig); end +fprintf('Complete data have been created.'); + +% Generation of missing array + +miss_array = cell(3*N,m); +vide = cell(1,2*N); l= 1; +while l <= m, + ev(1:N) = data(:,l); ev(N+1:3*N) = vide; + miss_array(:,l) = sample_bnet(bnet_miss, 'evidence', ev); + % apply simple rule of missingness + ev2 = cell2mat(miss_array(N+1:2*N, l)); + % verification that we have not a completly missing sample + ev2 = 3-ev2; + if prod(ev2)==1, + if v, fprintf(' - %d, one completly missing sample removed', l);end + else l=l+1; + end + if v, if mod(l,250)==0, fprintf('\n - %d',l); end, end +end +fprintf('\n'); + +% Generation of incomplete dataset + + %% TO REPLACE THE CELL ARRAY FOR ouput data + %% WITH A MATRIX WITH A SPECIAL CASE (size+1) + %% FOR MISSING DATA, SIMPLY REPLACE 1 by 0 +if 1, %% HERE + miss_array = bnt_to_mat(miss_array(N+1:2*N, :)); + miss_array = 2-miss_array; + data = bnt_to_mat(data); + comp_data = data; + data = data.*miss_array; + data = mat_to_bnt(data, 0); +else + comp_data = bnt_to_mat(miss_array(1:N,:),0); + data = bnt_to_mat(miss_array(2*N+1:3*N,:),0); + miss_array = bnt_to_mat(miss_array(N+1:2*N, :)); + miss_array = 2-miss_array; +end +fprintf('Incomplete dataset have been created.\n'); + +% Verification of the Rate of missing data + +if base_proba, + [XX, YY]=find(miss_array==0); + nbr_miss = length(YY); + taux = nbr_miss/N/m; + if v, fprintf('There is %2.2f percent of missing data\n', round(taux*10000)/100); end + + % Khi2 test between taux and base_proba for m*N + toto = m*N; + d = ((nbr_miss-base_proba*toto)^2)/(base_proba*toto) + (((toto-nbr_miss)-(1-base_proba)*toto)^2)/((1-base_proba)*toto); + if d>choice, + fprintf('THE DATASET DO NOT RESPECT %2.1f%% OF MISSING DATA (%2.1f%%, Khi2 : %2.1f > %2.1f)\n', round(base_proba*10000)/100, round(taux*10000)/100, d, choice); + notok = 1; + end +end diff --git a/sourcecodes/bnt-master/SLP/misc/gener_discrete_dist.m b/sourcecodes/bnt-master/SLP/misc/gener_discrete_dist.m new file mode 100644 index 00000000..4f131fcb --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/gener_discrete_dist.m @@ -0,0 +1,12 @@ +function beta = gener_discrete_dist(N, alpha) +% beta = gener_beta(N, masse) +% +% using the logistic function to create a random vector +% of discrete probability distribution on N values +% if masse~=1 then sum(beta)==masse and this is not a probability +% +% francois.olivier.c.h@gmail.com + +if nargin==1, alpha=1; end +gamma = rand(1,N); +beta = alpha*N*(exp(gamma)/sum(exp(gamma))); diff --git a/sourcecodes/bnt-master/SLP/misc/gener_empty_cache.m b/sourcecodes/bnt-master/SLP/misc/gener_empty_cache.m new file mode 100644 index 00000000..6e9669e5 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/gener_empty_cache.m @@ -0,0 +1,26 @@ +function cache = gener_empty_cache(N,L) +% cache = gener_empty_cache(number_of_nodes,lenght_of_cache) +% +% exemple for 2 nodes with cache of size 5 : +% +% cache = +% 5 b 0 0 0 --> 1st empty place 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 +% | | | | | +% | | | | |___> 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 +% +% +% +% designed by francois.olivier.c.h@gmail.com +% + +cache=zeros(L+1,N+3); +cache(1,1)=2; \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/hist_ic.m b/sourcecodes/bnt-master/SLP/misc/hist_ic.m new file mode 100644 index 00000000..b16980fd --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/hist_ic.m @@ -0,0 +1,268 @@ +function [n,edges,nbedges,xechan] = hist_ic(x,crit) + +%HIST_IC optimal Histogram based on IC information criterion +% +% [N,EDGES,NBEDGES,XECHAN] = HIST-IC(X,CRIT) +% bins the elements of X into an optimal number of bins according +% to a cost function based on Akaike's Criterion. +% +% +% CRIT = 1 | 2 | 3 (choose one of the 3 possible criterium) (default 3) +% 4 (returns the initial histogram instead of the optimal one) +% +% +% N = cell array containing the distribution of each column of X +% (or a vector if X is a column vector) +% EDGES = cell array containing the bin edges of each column of X +% (or a vector if X is a column vector) +% NBEDGES = vector containing the number of bin edges for each column of X +% (or a number if X is a column vector) +% XECHAN = discretized version of X +% +% Ref : O. Colot et al., Information Criteria and Abrupt Changes in +% Probability Laws, Signal Processing VII: Theory and Applications +% pp.1855-18858, September 1994 +% +% F. El-Matouat, O. Colot 2000 (first version) +% Revised 01-06-2001 by Ph. Leray - philippe.leray@univ-nantes.fr +% +% +% Things to do : +% * Call criteron by a name ('aic','xxx', ...) instead of a number +% + + +if nargin == 0 + error('Requires one or two arguments.') +end + +if nargin == 1 + crit = 3; +end; + +if min(size(x))==1, x = x(:); end + +if isstr(x) + error('Input argument must be numeric.') +end + +if isempty(x), + error('No elements to count') +end + + +[nb_l,nb_c]=size(x); + +% Outputs declaration +xechan=zeros(nb_l,nb_c); + +edges=cell(nb_c,1); +% Local variables +maxi = max(x); +mini = min(x); + +%% Erreur ? ancien code : +%% nb_clas_ini=2*round(sqrt(nb_l)-1); % article Fatima + +nb_clas_ini=round(2*sqrt(nb_l)-1); + +pas_ini=(maxi-mini)/nb_clas_ini; % initial step + +for j=1:nb_c, + + % optimal histogram for each column of X + histo_ini =hist(x(:,j),nb_clas_ini); % initial histogram + + if (crit~=4) + [hist_opt,pas_opt]=hist1_ic(histo_ini,nb_l,pas_ini(j),nb_clas_ini,crit); + else + fprintf('Histo initial\n'); + hist_opt=histo_ini; + pas_opt=ones(1,nb_clas_ini)*pas_ini(j); + end; + nbedges(j)=size(hist_opt,2); + edges{j}=mini(j)+cumsum(pas_opt(1:nbedges(j)-1)); %+1e-7; + [n{j} xechan(:,j)]=histc(x(:,j),[-inf edges{j} inf]); + n{j}=n{j}(1:end-1); +end + +if (nb_c==1) + n=n{1}; edges=edges{1}; +end + + +% ============================== subfunctions + +function [hist_opt,step_opt]=hist1_ic(histo,nb,step_ini,m,critere); + +%HIST1_IC optimal Histogram based on IC information criterion +% +% [HIST_OPT, STEP_OPT] = HIST1_IC(HISTO, NB, STEP_INI, NBSTEP_INI, CRIT) +% fusion of an 1D histogramme (HISTO) according to an IC criterion (CRIT) +% +% This function is mainly an internal function used by HIST_IC +% +% Ref : O. Colot et al., Information Criteria and Abrupt Changes in +% Probability Laws, Signal Processing VII: Theory and Applications +% pp.1855-18858, September 1994 +% +% F. El-Matouat, O. Colot 2000 (first version) +% Revised 11-06-2001 by Ph. Leray +% +% +% Things to do : +% * Call criteron by a name ('aic','xxx', ...) instead of a number +% + + +aic=[]; +aic2=[]; + +% Initialisation +histt=histo; +teta = histt/nb; +pas = step_ini*ones(1,m); + +% Calcul de l'ensemble des histogrammes optimaux + +for z=1:m + + % Calcul de AIC pour l'union entre hist(indice,u) et hist(indice,u+1) + aic2 = [aic2 cal_aic(nb,teta,pas,m+1-z,critere)]; + + if (z~=m) + % Calcul des couples de classes adjacentes + if critere==1 + penalite=(2*(m-z)-1)/nb; + elseif critere==2 + penalite=(m-z-1)*(1+log(nb))/nb; + else + penalite=(m-z)*(1+log(log(nb)))/nb; + end + aic=cla_adj(nb,teta,pas,step_ini,m-z+1,histt,penalite,aic); + % Recherche de la valeur min du crit�re pour les classes adjacentes + [min_aic classe]=min(aic(1:(m-z))); + + + % Fusion de hist(classe) et hist(classe+1) + nb_pas1=pas(classe)/step_ini; + nb_pas2=pas(classe+1)/step_ini; + + ess=round( nb_pas1*histt(classe)+nb_pas2*histt(classe+1) ); + teta(classe)=ess / nb; + histt(classe)=ess / (nb_pas1+nb_pas2); + pas(classe)=pas(classe)+pas(classe+1); + + + % Cr�ation du nouvel histogramme + itemp = setdiff(1:m+1-z,classe+1); + histt = histt(itemp); + pas = pas(itemp); + teta = teta(itemp); + end +end + +% Recherche du crit�re minimun AIC +[min_AIC fusion]=min(aic2(1:m)); + +% Initialisation de histo +histt=histo; +teta = histt/nb; +pas = step_ini*ones(1,m); + +% Calcul de l'histogramme optimal + +for z=1:fusion-1 + + % Calcul des couples de classes adjacentes + if critere==1 + penalite=(2*(m-1)-1)/nb; + elseif critere==2 + penalite=(m-2)*(1+log(nb))/nb; + else + penalite=(m-1)*(1+log(log(nb)))/nb; + end + + aic=cla_adj(nb,teta,pas,step_ini,m,histt,penalite,aic); + + + % Recherche de la valeur min du crit�re pour les classes adjacentes + [min_aic classe]=min(aic(1:m-1)); + + % Fusion de hist(indice,classe) et hist(indice,classe+1) + + nb_pas1=pas(classe)/step_ini; + nb_pas2=pas(classe+1)/step_ini; + + teta(classe)=(round(nb_pas1*histt(classe)+nb_pas2*histt(classe+1)))/nb; + histt(classe)=(nb_pas1*histt(classe)+nb_pas2*histt(classe+1))/(nb_pas1+nb_pas2); + pas(classe)=pas(classe)+pas(classe+1); + + % Cr�ation du nouvel histogramme + + itemp=setdiff(1:m,classe+1); + histt = histt(itemp); + pas = pas(itemp); + teta = teta(itemp); + %aic=zeros(1,m-1); + + m=m-1; +end +hist_opt=histt; +step_opt=pas; + + +%===================================================== +% Calcul du Critere pour l'ensemble des classes + +function akaike=cal_aic(size_ech,teta,pas,m,critere); + + +if critere==1 + a=(2*m-1)/size_ech; +elseif critere==2 + a=(m-1)*(1+log(size_ech))/size_ech; +else + a=m*(1+log(log(size_ech)))/size_ech; +end + +indu = find(teta); +akaike = a - 2*sum(teta(indu).*log(teta(indu)./pas(indu))); + + +%===================================================== +% Cla_adj.m +% aic=cla_adj(taille,indice,teta,pas,pas_ini,m,hist,penalite,aic) +% taille=nombre d'�l�ments dans chacune des hypotheses; +% indice=numero de la classe; +% Calcul du critere de Akaike pour l'histogramme totale avec +% fusion de deux classes adjacentes u et (u+1). + +function aic=cla_adj(size_ech,teta,pas,pas_ini,m,hist,penalite,aic); + + +for u=1:m-1 + + cumul=0; + + % This loop is faster than a sum of a vectorised computation ! + for x=1:m + if x~=u & x~=u+1 & teta(x)~=0 + cumul=cumul+teta(x)*log(teta(x)/pas(x)); + end + end + + + nb_pas1=pas(u)/pas_ini; + nb_pas2=pas(u+1)/pas_ini; + + b=( round(nb_pas1*hist(u)+nb_pas2*hist(u+1) ) ) / size_ech; + + if b~=0 + c=2*b*log( b / ( pas(u) + pas(u+1) ) ); + else + c=0; + end + + aic(u)=penalite-2*cumul-c; +end diff --git a/sourcecodes/bnt-master/SLP/misc/histc_ic.m b/sourcecodes/bnt-master/SLP/misc/histc_ic.m new file mode 100644 index 00000000..3e4a5af5 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/histc_ic.m @@ -0,0 +1,35 @@ +function [n,xechan] = histc_ic(x,edges) +%HISTC_IC Histogram count +% +% [N XECHAN] = HISTC_IC(X,EDGES), for vector X, counts the number of values +% in X that fall between the elements in the EDGES vector +% (EDGES = cell array returned by HIST_IC function) +% +% N is cell array containing these counts. +% (or a vector if X is a column vector) +% +% XECHAN = discretized version of X +% +% Example : +% X=randn(100,5); +% Xapp=X(1:50,:); +% Xtest=X(51:100,:); +% +% % bins are computed with Xapp data +% [n1,bornes,nbbornes,xechan]=hist_ic(Xapp); +% % histogram is computed on Xtest data, with previously obtained bins +% [n2,xtechan]=histc_ic(Xtest,bornes); +% +% 05-06-2001 by Ph. Leray - philippe.leray@univ-nantes.fr +% + +[nb_l,nb_c]=size(x); +% Outputs declaration +xechan=zeros(nb_l,nb_c); +if nb_c==1, edges={edges}; end +% edges=mat2cell(edges); modified by francois.olivier.c.h@gmail.com +for j=1:nb_c, + [n{j} xechan(:,j)]=histc(x(:,j),[-inf edges{j} inf]); + n{j}=n{j}(1:end-1); +end +if nb_c==1, n=n{1}; end diff --git a/sourcecodes/bnt-master/SLP/misc/inference.m b/sourcecodes/bnt-master/SLP/misc/inference.m new file mode 100644 index 00000000..3aef9eeb --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/inference.m @@ -0,0 +1,25 @@ +function [proba_post, engin2]= inference(bnet, data, node) +% Make bayesian inference on data +% [proba_post, engine]= inference(bnet, data, node) +% +% INPUTS : +% - bnet, the structure of the bayesian network gived by mk_bnet. +% - data(i,m), node i in case m. +% - node, the node we interrogating. +% +% OUTPUTS : +% - proba_post, the posteriors probabilities. +% - engine, the inference engine. +% +% francois.olivier.c.h@gmail.com + +engine=jtree_inf_engine(bnet); +[N L]=size(data); +proba_post=zeros(L,bnet.node_sizes(node)); +for i=1:L + evidence(1:N)=data(1:N,i); + evidence{node}=[]; + [engin2, ll]=enter_evidence(engine,evidence); + marg=marginal_nodes(engin2,node); + proba_post(i,:)=marg.T'; +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/isdag.m b/sourcecodes/bnt-master/SLP/misc/isdag.m new file mode 100644 index 00000000..8335d58c --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/isdag.m @@ -0,0 +1,9 @@ +function b = isdag(G) +% b = isdag(G) +% +b = sum(sum(G.*G')); % How many undirected arcs ? (x2) +b=~b & ~isempty(G); +if b + M = expm(full(G)) - eye(length(G)); M = (M>0); + b = isempty(find(sum(sum(eye(length(G)).*M)))); % is there no cycle ? +end diff --git a/sourcecodes/bnt-master/SLP/misc/ismemberclique.m b/sourcecodes/bnt-master/SLP/misc/ismemberclique.m new file mode 100644 index 00000000..5b1f8771 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/ismemberclique.m @@ -0,0 +1,14 @@ +function resu = ismemberclique(v,cliques) +% b = ismemberclique(v,cliques) +% + +finiclique = 0 ; resu=0 ; +cl=1; ncl=length(cliques) ; + +while (~finiclique) & (cl<=ncl); + if ismember(v,cliques{cl}) + resu=1; + finiclique=1 ; + end + cl=cl+1; +end diff --git a/sourcecodes/bnt-master/SLP/misc/knn.m b/sourcecodes/bnt-master/SLP/misc/knn.m new file mode 100644 index 00000000..43233de1 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/knn.m @@ -0,0 +1,57 @@ +function [ypred]=knn(xapp,yapp,valY,X,k) + +% +% knn implementation +% +% USE : [ypred]=knn(xapp,yapp,valY,X,k) +% +% Vincent Guigue 08/01/03 + +% check nargin + +if nargin<4 + error('too few argumemnts'); +elseif nargin<5 + k=3; +else + if mod(k,2)==0 + error('k must be odd'); + end +end + +if size(xapp,2)~=size(X,2) + error('dimension incompatibility'); +end + + +ndim = size(xapp,2); +nptxapp = size(xapp,1); +nptX = size(X,1); + +% distance de X a xapp : +mat1 = repmat(xapp, nptX,1); +%mat21 = reshape(X',1,nptX*ndim) +mat22 = repmat(X,1,nptxapp)'; +mat2 = reshape(mat22 ,ndim, nptxapp*nptX)'; +distance = mat1 - mat2 ; + +distance = sum(distance.^2,2); +distance = reshape(distance,nptxapp,nptX); +[val kppv] = sort(distance,1); + +% bilan sur les k premieres lignes +kppv = reshape(kppv(1:k,:),k*nptX,1); +Ykppv = yapp(kppv,1); +Ykppv = reshape(Ykppv,k,nptX); + +% trouver le plus de reponses identique par colonne + +vote = []; +for i=1:nptX + for j=1:length(valY) + vote(j,i)=size(find(Ykppv(:,i)==valY(j)),1); + end +end + +[val ind]=max(vote,[],1); +ypred = valY(ind); \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/mat_to_bnt.m b/sourcecodes/bnt-master/SLP/misc/mat_to_bnt.m new file mode 100644 index 00000000..008cabc8 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/mat_to_bnt.m @@ -0,0 +1,58 @@ +function [res] = mat_to_bnt(mat,misv) +% MAT_TO_BNT Convert a matrix to a cell array +% D = mat_to_bnt(data,misv) +% +% Input : +% data(i,m) is the node i in the case m, +% misv is the way you choose to encode missing data +% in the original matrix (-9999 by default) +% +% Output : +% D = cell array containing data(i,m) is the data is OK +% or [] is the data is missing +% +% +% V1.2 : 18 feb 2003 (Ph. Leray - philippe.leray@univ-nantes.fr) +% +% >> m=rand(2,4) +% +% m = +% +% 0.9525 0.4693 0.3907 0.1496 +% 0.9274 0.3157 0.1346 0.9383 +% +% >> m(2,2)=-9 +% +% m = +% +% 0.9525 0.4693 0.3907 0.1496 +% 0.9274 -9.0000 0.1346 0.9383 +% +% >> mat_to_bnt(m,-9) +% +% ans = +% +% [0.9525] [0.4693] [0.3907] [0.1496] +% [0.9274] [] [0.1346] [0.9383] +% + +if nargin <1 + error('Requires at least 1 argument.') +end + +if nargin == 1 + misv=-9999; + end; + +taille=size(mat); +long=taille(1); +larg=taille(2); +for i=1:long + for j=1:larg + res{i,j}=mat(i,j); + if(mat(i,j)==misv) + res{i,j}=[]; + end + end +end + diff --git a/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_dag_topo.m b/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_dag_topo.m new file mode 100644 index 00000000..73ba2439 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_dag_topo.m @@ -0,0 +1,196 @@ +function [new_nbrs, new_ops, new_nodes, new_topos] = mk_nbrs_of_dag_topo(G0) +% MK_NBRS_OF_DAG_TOPO Make all DAGs that differ from G0 by a single edge deletion, addition or reversal +% [new_nbrs, new_ops, new_nodes, new_topos] = mk_nbrs_of_dag_topo(G0) +% +% new_nbrs{i} is the i'th neighbor of G0. +% new_ops{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor. +% new_nodes(i,1:2) are the head and tail of the operated-on arc. +% new_topos are topological orders of the neighbours. +% +% We implement the fast acyclicity check described by P. Giudici and R. Castelo, +% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001. +% +% Written by Qian Diao <qian.diao@intel.com> on 19 Nov 01 +% Reference are ..\BNT\graph\mk_nbrs_of_dag.m, ..\BNT\learning\learn_struct_mcmc.m and +% ..\BNT\graph\topological_sort.m +% Copyright Intel 2001 +% +new_nbrs = {}; +new_ops = {}; +new_nodes = []; +new_topos = {}; +cs = {}; + +n = length(G0); +indeg = zeros(1,n); +zero_indeg = []; % a stack of nodes with no parents +for i=1:n + indeg(i) = length(parents(G0,i)); + cs{i} = children(G0, i); + if indeg(i)==0 + zero_indeg = [i zero_indeg]; + end +end + +dag = G0; +[nbrs, ops, nodes] = mk_nbrs_of_digraph(dag); +A = init_ancestor_matrix(dag); +%assert(acyclic(new_dag)); + +d1 = 1; +for d = 1:length(ops) + i = nodes(d, 1); j = nodes(d, 2); + legal = 0; + switch ops{d} + case 'add', + if A(i,j)==0 + legal = 1; + end + case 'del', + legal = 1; + case 'rev', + ps = mysetdiff(parents(dag, j), i); + % if any(A(ps,i)) then there is a path i -> parent of j -> j + % so reversing i->j would create a cycle + legal = ~any(A(ps, i)); + end + + if legal + tmp = nbrs(:,:,d); + new_nbrs{d1} = tmp; + new_ops{d1} = ops{d}; + new_nodes(d1,1:2) = nodes(d,1:2); + + % obtain the topological orders of neighbour dags + zero_indeg_nbr = []; + indeg_nbr = []; + cs_nbr = []; + + switch ops{d} + case 'add' % i is a new parent of j + zero_indeg_nbr = zero_indeg; + indeg_nbr = indeg; + if ~isempty(find(zero_indeg == j)) % j is not a root anymore + zero_indeg_nbr = mysetdiff(zero_indeg, j); + end + indeg_nbr(j) = indeg(j)+1; + + t_nbr=1; + order_nbr = zeros(1,n); + while ~isempty(zero_indeg_nbr) + v_nbr = zero_indeg_nbr(1); % pop v + zero_indeg_nbr = zero_indeg_nbr(2:end); + order_nbr(t_nbr) = v_nbr; + t_nbr = t_nbr + 1; + if v_nbr == i % j is a new child of i + cs_nbr = sort([j cs{i}]); + else + cs_nbr = cs{v_nbr}; + end + for k = 1:length(cs_nbr) + c_nbr = cs_nbr(k); + indeg_nbr(c_nbr) = indeg_nbr(c_nbr) - 1; + if indeg_nbr(c_nbr) == 0 + zero_indeg_nbr = [c_nbr zero_indeg_nbr]; % push c + end + end + end + + case 'del' % i is not a parent of j anymore + zero_indeg_nbr = zero_indeg; + indeg_nbr = indeg; + if length(parents(tmp, j))==0 + zero_indeg_nbr = -sort(-[zero_indeg, j]); % descending order + end + indeg_nbr(j) = indeg(j) - 1; + + t_nbr=1; + order_nbr = zeros(1,n); + while ~isempty(zero_indeg_nbr) + v_nbr = zero_indeg_nbr(1); % pop v + zero_indeg_nbr = zero_indeg_nbr(2:end); + order_nbr(t_nbr) = v_nbr; + t_nbr = t_nbr + 1; + if v_nbr == i % j is not a child of i anymore + cs_nbr = mysetdiff(cs{i}, j); + else + cs_nbr = cs{v_nbr}; + end + for k = 1:length(cs_nbr) + c_nbr = cs_nbr(k); + indeg_nbr(c_nbr) = indeg_nbr(c_nbr) - 1; + if indeg_nbr(c_nbr) == 0 + zero_indeg_nbr = [c_nbr zero_indeg_nbr]; % push c + end + end + end + + case 'rev' %i is a new child of j and j is a new parent of i + zero_indeg_nbr = zero_indeg; + indeg_nbr = indeg; + if ~isempty(find(zero_indeg == i)) + zero_indeg_nbr = mysetdiff(zero_indeg_nbr, i); + end + if length(parents(tmp, j))==0 + zero_indeg_nbr = -sort(-[zero_indeg_nbr, j]); % decending order + end + indeg_nbr(i) = indeg(i)+1; + indeg_nbr(j) = indeg(j)-1; + + t_nbr=1; + order_nbr = zeros(1,n); + while ~isempty(zero_indeg_nbr) + v_nbr = zero_indeg_nbr(1); % pop v + zero_indeg_nbr = zero_indeg_nbr(2:end); + order_nbr(t_nbr) = v_nbr; + t_nbr = t_nbr + 1; + cs_nbr = cs{v_nbr}; + if v_nbr == i % j is not a child of i anymore + cs_nbr = mysetdiff(cs{i}, j); + end + if v_nbr == j % i is a new child of j + cs_nbr = sort([i cs{j}]); + end + for k = 1:length(cs_nbr) + c_nbr = cs_nbr(k); + indeg_nbr(c_nbr) = indeg_nbr(c_nbr) - 1; + if indeg_nbr(c_nbr) == 0 + zero_indeg_nbr = [c_nbr zero_indeg_nbr]; % push c + end + end + end + end + + new_topos{d1} = order_nbr; + d1 = d1+1; + end +end + +clear nbrs ops nodes; + + + +%%%%%%%%% +function A = update_row(A, j, dag) +% We compute row j of A +A(j, :) = 0; +ps = parents(dag, j); +if ~isempty(ps) + A(j, ps) = 1; +end +for k=ps(:)' + anck = find(A(k,:)); + if ~isempty(anck) + A(j, anck) = 1; + end +end + +%%%%%%%% +function A = init_ancestor_matrix(dag) +order = topological_sort(dag); +A = zeros(length(dag)); +for j=order(:)' + A = update_row(A, j, dag); +end + + \ No newline at end of file diff --git a/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_pdag_add.m b/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_pdag_add.m new file mode 100644 index 00000000..f495702b --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_pdag_add.m @@ -0,0 +1,169 @@ +function [PDAGs, nodes] = mk_nbrs_of_pdag_add(cpdag,engine) +% MK_NBRS_OF_PDAG_ADD Make the superior inclusion boundary of CPDAG. +% [PDAGs, nodes] = mk_nbrs_of_pdag_add(CPDAG) +% +% PDAGs{i} is the i'th neighbor of CPDAG0 generated by INSERT(X,Y,T) with +% nodes{i,1:2}=[X Y] +% nodes{i,3}=T +% +% See D.M. Chickering 2002 : "Optimal Structure Identification with Greedy Search". +% +% philippe.leray@univ-nantes.fr, francois.olivier.c.h@gmail.com +% 24 july 2003 + + +compteur=0 ; +N=length(cpdag); +G=pdag_to_dag(cpdag); +if nargin==1, + bnet_tmp=mk_bnet(pdag_to_dag(cpdag),2*ones(N,1)); + bnet_tmp=mk_bnet(G,2*ones(N,1)); + engine_tmp=struct(jtree_inf_engine(bnet_tmp)); + clear bnet_tmp +else + engine_tmp=struct(engine); +end +cliques=engine_tmp.cliques; +nbcliques=length(cliques); +clear engine_tmp; + +% find in the PDAG all the X Y not connected +[LX LY]=find((cpdag|cpdag')+eye(N)==0); +nlinks=length(LX); % here is a bug when nlinks is zeros, for i=1:nlinks fail 31-37 added by hanbin +if nlinks==0 + [i,j]=find(cpdag);k=unidrnd(size(i,1),1); + cpdag(i(k),j(k))=0; + [LX LY]=find((cpdag|cpdag')+eye(N)==0); + nlinks=length(LX); +end + +for i=1:nlinks + X=LX(i); + Y=LY(i); + % Neighbors of Y + NY = myintersect(find(cpdag(:,Y)), find(cpdag(Y,:))); + % Adjacents of X + AX = myunion(find(cpdag(:,X)), find(cpdag(X,:))); + % Neighbors of Y adjacent to X + NAYX = myintersect(NY,AX); + % Neighbors of Y NOT adjacent to X + NNAYX = mysetdiff(NY, NAYX); + + % this function recursively "walks" (dfs) in the graph representation of NNA powerset + liste=NNAYX; + if ~isempty(liste) + premier=liste(1); + dernier=liste(end); + end; + + current_set=[]; + fini=0; + evite2 = 0 ; + + while ~fini + isclique=0; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 1 + % is [NAYX current_set] a clique ? + if isempty(current_set) + NAYXT=NAYX; + elseif isempty(NAYX) + NAYXT=current_set; + else + NAYXT = union(NAYX,current_set); + end + % if isempty(NAYXT),isclique=1;end + + isclique = isempty(NAYXT) | ismemberclique(NAYXT,cliques) ; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 1 + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 2 + % is there exist a partially directed path Y...X in PDAG \ NAYXT ? + if isclique + if evite2 + test2=1 ; + else + %%%% calcul test 2 + L2 = setdiff(1:N,NAYXT); + test2=~partialconnected(cpdag(L2,L2),find(L2==Y),find(L2==X)); + + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 2 + + if test2 + evitet2=1; + + % avoid testing INSERT(X,Y,[]) and INSERT(Y,X,[]) if pa(Y)=pa(X) + test0=1 ; + if (X>Y) + if length(current_set)==0 + PaX=setdiff(find(cpdag(X,:)),find(cpdag(:,X))); + PaY=setdiff(find(cpdag(Y,:)),find(cpdag(:,Y))); + test0=(length(setdiff(PaX,PaY))~=0); + end + end + if test0 + %fprintf(' INSERT(%d,%d,',X,Y); fprintf('%d',current_set); + %fprintf(')\n'); + compteur=compteur+1; + nodes{compteur,1}=X; + nodes{compteur,2}=Y; + nodes{compteur,3}=current_set; + ptmp=cpdag; + ptmp(X,Y)=1; + ptmp(current_set,Y)=1; + ptmp(Y,current_set)=0; + PDAGs{compteur}=ptmp; + end + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next Set to test ... + % what is the net set in the powerset ? + if length(liste)==0 + fini=1 ; next_set=[]; + elseif length(current_set)==0 + next_set=[premier]; % first node after the root [] + else + actuel=current_set(end); + if actuel==dernier + if length(current_set)==1 + fini=1; % no more node ... + else + ancien=current_set(end-1); % new "branch" + next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)]; + end + else % new node in the "branch" + if ~isclique + if length(current_set)==1 + fini=1; % no more node ... + else + ancien=current_set(end-1); % new "branch" + next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)]; + evite2=0; + end + else + next_set=[current_set liste(find(liste==actuel)+1)]; + end + end + end + current_set=next_set; + end +end + + +%%%%%%%%% +function resu = ismemberclique(v,cliques) +finiclique = 0 ; resu=0 ; +cl=1; ncl=length(cliques) ; +while (~finiclique) & (cl<=ncl); + if ismember(v,cliques{cl}) + resu=1; + finiclique=1 ; + end + cl=cl+1; +end + +%%%%%%%%% +function resu = partialconnected(G,Y,X) +tmp=expm(G); +resu = (tmp(Y,X)~=0); diff --git a/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_pdag_del.m b/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_pdag_del.m new file mode 100644 index 00000000..7f0fc35a --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/mk_nbrs_of_pdag_del.m @@ -0,0 +1,175 @@ +function [PDAGs, nodes] = mk_nbrs_of_pdag_del(cpdag,engine) +% MK_NBRS_OF_PDAG_ADD Make the inferior inclusion boundary of CPDAG. +% [PDAGs, nodes] = mk_nbrs_of_pdag_del(CPDAG) +% +% PDAGs{i} is the i'th neighbor of CPDAG0 generated by DELETE(X,Y,H) with +% nodes{i,1:2}=[X Y] +% nodes{i,3}=H +% +% See D.M. Chickering 2002 : "Optimal Structure Identification with Greedy Search". +% +% philippe.leray@univ-nantes.fr +% 25 july 2003 + +compteur=0; +N=length(cpdag); +if nargin==1, + bnet_tmp=mk_bnet(pdag_to_dag(cpdag),2*ones(N,1)); + engine_tmp=struct(jtree_inf_engine(bnet_tmp)); + clear bnet_tmp +else + engine_tmp=struct(engine); +end +cliques=engine_tmp.cliques; +nbcliques=length(cliques); +clear engine_tmp; + +verbose=0; + +if verbose + fprintf('---- Les cliques (maximales) du DAG\n'); + for i=1:nbcliques, + disp(cliques{i}); + end + fprintf('---- Variables \n'); +end + + +% find in the PDAG all the X Y connected +[LX LY]=find(cpdag==1); +nlinks=length(LX); + +for i=1:nlinks + X=LX(i); + Y=LY(i); + % Neighbors of Y + NY = myintersect(find(cpdag(:,Y)), find(cpdag(Y,:))); + % Adjacents of X + AX = myunion(find(cpdag(:,X)), find(cpdag(X,:))); + % Neighbors of Y adjacent to X + NAYX = myintersect(NY,AX); + + % this function recursively "walks" (dfs) in the graph representation of NA powerset + if verbose + X, Y, NY, AX, NAYX + fprintf('---- NA Powerset\n\n'); + end + + liste=NAYX; + if ~isempty(liste) + premier=liste(1); + dernier=liste(end); + end + + current_set=[]; + fini=0; + evite2 = 0 ; + + while ~fini + if verbose + % sert uniquement � l'affichage pour le prog de test ... + if isempty(current_set) + fprintf('\t H = []\n\n'); + else + fprintf('\t H = '); disp(current_set); + end + end + + isclique=0; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 1 + % is [NAYX \ current_set] a clique ? + if isempty(current_set) + NAYXH=NAYX; + elseif isempty(NAYX) + NAYXH=current_set; + else + NAYXH = mysetdiff(NAYX,current_set); + end + + if isempty(NAYXH) + if verbose + fprintf('\t\t NA(Y,X) \\ T = \t[]\n\n'); + end + + elseif verbose + fprintf('\t\t NA(Y,X) \\ T ='); disp(NAYXH); + end + + isclique = isempty(NAYXH) | ismemberclique(NAYXH,cliques) ; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 1 + if isclique + % avoid testing INSERT(X,Y,H) and INSERT(Y,X,H) if X--Y + test0=1 ; + if (X>Y) + test0=~cpdag(Y,X); + if verbose&~test0 + fprintf(' (inutile de tester %d %d ',X,Y); fprintf('%d',current_set); + fprintf(')\n'); + end + end + if test0 + if verbose + fprintf(' ==DELETE(%d,%d,',X,Y); fprintf('%d',current_set); + fprintf(')\n'); + end + compteur=compteur+1; + nodes{compteur,1}=X; + nodes{compteur,2}=Y; + nodes{compteur,3}=current_set; + ptmp=cpdag; + ptmp(X,Y)=0; + ptmp(Y,X)=0; + ptmp(current_set,Y)=0; + ptmp(Y,current_set)=1; + ptmp(current_set(find(cpdag(X,current_set)==1)),X)=0; + PDAGs{compteur}=ptmp; + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next Set to test ... + % what is the net set in the powerset ? + if length(liste)==0 + fini=1 ; next_set=[]; + elseif length(current_set)==0 + next_set=[premier]; % first node after the root [] + else + actuel=current_set(end); + if actuel==dernier + if length(current_set)==1 + fini=1; % no more node ... + else + ancien=current_set(end-1); % new "branch" + next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)]; + end + else % new node in the "branch" + if ~isclique + if length(current_set)==1 + fini=1; % no more node ... + else + ancien=current_set(end-1); % new "branch" + next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)]; + evite2=0; + end + else + next_set=[current_set liste(find(liste==actuel)+1)]; + end + end + end + current_set=next_set; + end +end + + +%%%%%%% +function resu = ismemberclique(v,cliques) + +finiclique = 0 ; resu=0 ; +cl=1; ncl=length(cliques) ; + +while (~finiclique) & (cl<=ncl); + if ismember(v,cliques{cl}) + resu=1; + finiclique=1 ; + end + cl=cl+1; +end diff --git a/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.c b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.c new file mode 100644 index 00000000..3d347f65 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.c @@ -0,0 +1,137 @@ +/* multiply_one_marginals.c */ +/* the first operand can be a joint marginals of nodes set,*/ +/* but the second operand nust be a single node's marginal.*/ +/* and the result joint marginal has domain like [prhs[0].domain, prhs[2].domain]*/ +/* i.e. cat the second domain at the end of the first domain*/ +/* the third operands will be the eff_node_sizes */ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + mxArray *ptemp, *ptemp1; + double *xdom, *ydom, *xp, *yp, *zp, *pr, *pSizes; + int N, xElements, numElements; + int i, j, nxdom, ndim, ydim; + int *xdim, *s, *sx, *sy, *cpsx, *cpsy; + int *subs, *s1, *cpsx2, *cpsy2; + int rdims[2]; + const char *field_names[] = {"domain", "T", "mu", "Sigma"}; + + if (nrhs != 3) + mexErrMsgTxt("Incorrect number of inputs."); + + if (nlhs > 1) + mexErrMsgTxt("Too many output arguments."); + + if(mxIsEmpty(prhs[0])){ + plhs[0] = mxDuplicateArray(prhs[1]); + return; + } + + N = mxGetNumberOfElements(prhs[2]); + pSizes = mxGetPr(prhs[2]); + + ptemp = mxGetField(prhs[0], 0, "domain"); + nxdom = mxGetNumberOfElements(ptemp); + xdom = mxGetPr(ptemp); + ptemp = mxGetField(prhs[1], 0, "domain"); + ydom = mxGetPr(ptemp); + ndim = nxdom + 1; + + rdims[0] = 1; + rdims[1] = 1; + plhs[0] = mxCreateStructArray(2, rdims, 4, field_names); + ptemp = mxCreateDoubleMatrix(1, ndim, mxREAL); + mxSetField(plhs[0], 0, "domain", ptemp); + pr = mxGetPr(ptemp); + for(i=0; i<nxdom; i++){ + pr[i] = xdom[i]; + } + pr[ndim-1] = *ydom; + + xdim = (int *)malloc(sizeof(int)*nxdom); + for(i=0; i<nxdom; i++){ + xdim[i] = (int)pSizes[(int)xdom[i]-1]; + } + ydim = (int)pSizes[(int)*ydom - 1]; + + ptemp = mxGetField(prhs[1], 0, "T"); + yp = mxGetPr(ptemp); + ptemp = mxGetField(prhs[0], 0, "T"); + xp = mxGetPr(ptemp); + xElements = mxGetNumberOfElements(ptemp); + if(ydim == 1){ + ptemp1 = mxDuplicateArray(ptemp); + mxSetField(plhs[0], 0, "T", ptemp1); + free(xdim); + return; + } + numElements = xElements * ydim; + + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + s = (int *)malloc(sizeof(int)*ndim); + s1 = (int *)malloc(sizeof(int)*ndim); + *(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1; + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsx2 = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<nxdom; i++){ + subs[i] = 0; + sx[i] = xdim[i]; + sy[i] = 1; + s[i] = sx[i]; + s1[i] = s[i] - 1; + } + subs[ndim-1] = 0; + sx[ndim-1] = 1; + sy[ndim-1] = ydim; + s[ndim-1] = ydim; + s1[ndim-1] = s[ndim-1] - 1; + + for(i=0; i<ndim-1; i++){ + cpsx[i+1] = cpsx[i]*sx[i]--; + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsx2[i] = cpsx[i]*sx[i]; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]); + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + ptemp = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + mxSetField(plhs[0], 0, "T", ptemp); + zp = mxGetPr(ptemp); + + for(j=0; j<numElements; j++){ + *zp++ = *xp * *yp; + for(i=0; i<ndim; i++){ + if(subs[i] == s1[i]){ + subs[i] = 0; + if(sx[i]) + xp -= cpsx2[i]; + if(sy[i]) + yp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sx[i]) + xp += cpsx[i]; + if(sy[i]) + yp += cpsy[i]; + break; + } + } + } + free(xdim); + free(sx); + free(sy); + free(s); + free(s1); + free(cpsx); + free(cpsy); + free(subs); + free(cpsx2); + free(cpsy2); +} diff --git a/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.dll b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.dll new file mode 100644 index 00000000..7277cda4 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.dll Binary files differdiff --git a/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.mexa64 b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.mexa64 new file mode 100644 index 00000000..b07db577 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.mexa64 Binary files differdiff --git a/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.mexglx b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.mexglx new file mode 100644 index 00000000..abb17b4d --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/multiply_one_marginal.mexglx Binary files differdiff --git a/sourcecodes/bnt-master/SLP/misc/pdag_to_all_dags.m b/sourcecodes/bnt-master/SLP/misc/pdag_to_all_dags.m new file mode 100644 index 00000000..addfb4b4 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/pdag_to_all_dags.m @@ -0,0 +1,123 @@ +function [n_dags,dag_list] = pdag_to_all_dags( pdag ) + +% +% [n_dags,dag_list] = pdag_to_all_dags( pdag) +% +% generates a cell array of ALL Markov-equivalent DAGs +% corresponding to a partially directed acyclic graph (PDAG). +% +% Input: PDAG (PDAG does NOT have to be complete) +% Format of pdag matrix: +% Edge with known direction a->b represented as pdag(a,b)=-1 pdag(b,a)=0 +% Edge with unknown direction a-b represented as pdag(a,b)=1 pdag(b,a)=1 +% +% Output: Number of DAGs generated and +% Cell array of all permissible extensions of PDAG +% +% Sample Use: +% % Use output of PC algorithm +% dag = mk_rnd_dag(4); % create random DAG +% % Generate pdag through PC algorithm +% pdag = learn_struct_pdag_pc('dsep', length(dag), 3, dag) +% [n_dags,dag_list] = pdag_to_all_dags(pdag); +% n_dags +% +% If you want to generate all DAGs that are Markov equivalent to an +% input DAG (not a pattern), then use function Markov_equivalent_dags(dag) +% instead which calls this function. +% +% ======================================================================= +% Algorithm to generate ALL DAGs (pdag_to_all_dags): +% +% 0) Initialize an empty list of DAGs. +% +% 1) Complete current PDAG as far as possible using Rules R1-R4. +% +% 2) Select an unoriented edge X-Y. +% +% a) If none left: +% Done. Add DAG=abs(PDAG) to list of output DAGs. Return. +% +% b) Otherwise: +% Select an unoriented edge X-Y. +% Create PDAG1 with X->Y. +% Create PDAG2 with Y->X. +% Recursion: For EACH PDAG (PDAG1/2): Go to Step 1. +% +% +% This algorithm is a slight modification of the algorithm by Meek (1995) +% which generates a single DAG extension of a PDAG - here we just add +% recursion to consider both possible orientations for each considered edge. +% +% For the original algorithm by Meek, see +% C. Meek, "Causal inference and causal explanation with background +% knowledge", UAI 1995, Section 3.1.1, "Phase III" algorithm. +% +% Thanks to Daniel Eaton for extensive testing and bug reports. +% +% Imme Ebert-Uphoff (ebert@tree.com), 2007 +% ======================================================================= + + % MAIN + dag_list={}; % init empty list of DAGs + + % Complete pdag as far as possible using Rules R1-R4 of Meek (1995) + cpdag = complete_pattern(pdag); + + % Start recursion + dag_list = recurse_unoriented_edge(cpdag,dag_list); + + % return # of DAGs along with dag_list + n_dags = length(dag_list); + if (n_dags == 0) % no DAGs generated + fprintf('\nPDAG does not have any permissible extension!\n'); + end +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% RECURSE_UNORIENTED_EDGE % +% implements Step 2 of the pdag_to_all_dags algorithm. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function updated_list = recurse_unoriented_edge(cpdag, dag_list) + + % input must be a COMPLETE pdag + + [A,B] = find(cpdag==1); % find all undirected edges + + updated_list = dag_list; + + if isempty(A) % if no undirected edges left + % End of recursion reached. + % Convert all (-1) values to (1) to yield standard DAG, add DAG to list. + updated_list{end+1} = abs(cpdag); + + else + a = A(1); b = B(1); % choose first unoriented edge + % (any unoriented edge could be used here) + + % choose two different directions for edge and complete BOTH ! + % PDAG1: contains a -> b + %fprintf('PDAG1: %d -> %d\n',a,b); + pdag1 = cpdag; + pdag1(a,b) = -1; pdag1(b,a) = 0; + % complete as far as possible using rules R1-R4: + cpdag1 = complete_pattern(pdag1); + % Continue recursion on another unoriented edge + updated_list = recurse_unoriented_edge(cpdag1,updated_list); + + % PDAG1: contains b -> a + %fprintf('PDAG2: %d -> %d\n',b,a); + pdag2 = cpdag; + pdag2(a,b) = 0; pdag2(b,a) = -1; + % complete as far as possible using rules R1-R4: + cpdag2 = complete_pattern(pdag2); + % Continue recursion on another unoriented edge + updated_list = recurse_unoriented_edge(cpdag2, updated_list); + + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + diff --git a/sourcecodes/bnt-master/SLP/misc/pdag_to_dag.m b/sourcecodes/bnt-master/SLP/misc/pdag_to_dag.m new file mode 100644 index 00000000..00a73aac --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/pdag_to_dag.m @@ -0,0 +1,87 @@ +function G2 = pdag_to_dag(pdags) +% (also works with a cell array of pdags, returning a cell array of dags) +% dag = pdag_to_dag(pdag) +% +% cf Dor and Tarsi (1992) : +% A simple algorithm to construct a consistent extention of a partially oriented graph. +% +% francois.olivier.c.h@gmail.com + + +if ~iscell(pdags) + pdag=cell(1,1); + pdag{1}=pdags; +else + pdag=pdags; +end + +for da=1:length(pdag) + %fprintf('%d ',da) + G=pdag{da}; + G2 = G; A = G; + N = size(G,1); + empty_loop = 0; + + while ~isempty(find(A)) + [x x_y_undirected] = select_vertex(A); + + if x==0 + fprintf('pdag_to_dag error : This pdag does not admit any extension.\n'); + G2=[]; + break + end + G2(x,x_y_undirected) = 0; G2(x_y_undirected,x) = 1; + + A(x,:) = 0; + A(:,x) = 0; + end + dags{da}=G2; +end +if ~iscell(pdags) + G2=dags{1}; +else + G2=dags; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [sol, x_y] = select_vertex(G) +N = size(G,1); +sol = 0; +x = 0; +fini=0 ; +while ~fini + + x = x+1; + if x>N + fini=1; + else + beforex=find(G(:,x)); + afterx=find(G(x,:)); + + if ~(isempty(beforex)&isempty(afterx)) + x_y = myintersect(beforex,afterx); + beforex=mysetdiff(beforex,x_y); + afterx=mysetdiff(afterx,x_y); + + if isempty(afterx) % x is a sink + Ax= myunion(x_y,beforex); + for y=x_y + % Adjacents of y + Ay = myunion(find(G(:,y)), find(G(y,:))); + Ay = myunion(Ay,y); + if isempty(setdiff(Ax,Ay)) + fini=fini+1; + else + break; + end + end + if fini==length(x_y) + sol=x; fini=1; + else + fini=0; + end + end + end + end +end % while diff --git a/sourcecodes/bnt-master/SLP/misc/pdag_unsigned_to_signed.m b/sourcecodes/bnt-master/SLP/misc/pdag_unsigned_to_signed.m new file mode 100644 index 00000000..37b414f1 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/pdag_unsigned_to_signed.m @@ -0,0 +1,22 @@ +function signed_pdag = pdag_unsigned_to_signed(pdag) + +% Input: pdag with '1' in one place for every directed edge and +% '1' in two places for every undirected edge +% +% Output: pdag with '-1' in one place for every directed edge +% '1' in two places for every undirected edge +% +% This function is used by Markov_equivalent_dags(dag) to convert +% output of SLP function dag_to_cpdag(dag) +% to format required as input for pdag_to_all_dags(pdag). +% +% I'm sure there's a prettier way to code this! +% +% Imme Ebert-Uphoff (ebert@tree.com), 2007 +% + + undirected = ( (pdag+pdag')/2 == 1); % extract undirected eges + directed = pdag - undirected; % extract directed edges + signed_pdag = undirected - directed; % 1s for undirected, (-1)s for directed + + diff --git a/sourcecodes/bnt-master/SLP/misc/subsets1.m b/sourcecodes/bnt-master/SLP/misc/subsets1.m new file mode 100644 index 00000000..ca2bf46c --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/subsets1.m @@ -0,0 +1,45 @@ +function sub_s=subsets1(s,k) +% SUBSETS1 creates sub-sets of a specific from a given set +% SS = subsets1(S, k) +% +% S is the given set +% k is the required sub-sets size +% +% Example: +% +% >> ss=subsets1([1:4],3); +% >> ss{:} +% ans = +% 1 2 3 +% ans = +% 1 2 4 +% ans = +% 1 3 4 +% ans = +% 2 3 4 +% +% Written by Raanan Yehezkel, 2004 + +if k<0 % special case + error('subset size must be positive'); +elseif k==0 % special case + sub_s={[]}; +else + l=length(s); + ss={}; + if l>=k + if k==1 % Exit condition + for I=1:l + ss{I}=s(I); + end + else + for I=1:l + ss1=subsets1(s([(I+1):l]),k-1); + for J=1:length(ss1) + ss{end+1}=[s(I),ss1{J}]; + end + end + end + end + sub_s=ss; +end |