diff options
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 |