diff options
Diffstat (limited to 'sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m')
| -rw-r--r-- | sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m | 725 |
1 files changed, 725 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m b/sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m new file mode 100644 index 00000000..bddce565 --- /dev/null +++ b/sourcecodes/bnt-master/SLP/learning/learn_struct_bnpc.m @@ -0,0 +1,725 @@ +function [Phase_3, Phase_2, Phase_1, UPhase_3] = learn_struct_bnpc(data,node_sizes,epsilon,star) +% G = learn_struct_bnpc(Data,node_sizes,epsilon,star) +% +% Data(i,m) is node i in case m. +% node_sizes and epsilon are optionnals. +% star = 0 to use try_to_separate_B instead of try_to_separate_B_star +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. +% +% Things to do : rewrite function orient_edges ! +% ! sometimes it causes crashes ! +% +% V0.91 : 18 sept 2003 (olivier.francois@insa-rouen.fr) + +verbose=1; +%if nargin < 5, mwst=0; end +if nargin < 4, star=1; end +if nargin < 3, epsilon=0.05; end +if nargin < 2, node_sizes=max(data'); end + +if verbose + fprintf('================== phase I : \n'); +end +tmp1=cputime; +[Phase_1 II JJ score_mat score_mat2] = phaseI(data, node_sizes, epsilon); +tmp1=cputime-tmp1; + +if verbose + fprintf('Execution time : %2.5f\n',tmp1); + fprintf('\n================== phase II : \n'); +end +tmp1=cputime; +Phase_2 = phaseII(Phase_1, data, node_sizes, epsilon, II, JJ, score_mat); +tmp1=cputime-tmp1; + +if verbose + fprintf('Execution time : %2.5f\n',tmp1); + fprintf('\n================== phase III : \n'); +end +tmp1=cputime; +[Phase_3 UPhase_3] = phaseIII(Phase_2, data, node_sizes, epsilon, score_mat2, star); +tmp1=cputime-tmp1; +if verbose + fprintf('Execution time : %2.5f\n',tmp1); +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [G, II, JJ, score_mat, sc2] = phaseI(data,node_sizes,alpha) +% [G, II , JJ, score_mat, s2] = phaseI(data,node_sizes,epsilon) +% +% G is an acyclic graph +% [II JJ] is the list of important edges not processed in phase I (for phase II) +% score_mat is the mutual information score matrix +% +% data(i,m) is node i in case m. +% alpha is the significant level for CI tests ( default=0.05 ). +% node_sizes is the vector of sizes ( default=max(data') ). +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +% 0. +if nargin < 3, alpha=0.05; end +if nargin < 2, node_sizes=max(data'); end +[N m] = size(data); +score_mat = zeros(N); +edges=0; + +% 1. +G = zeros(N); +L=[]; + +% 2. Use of Chi2 instead of MI ... allow using a confidence level alpha instead of an arbitrary epsilon +for i=1:(N-1) + for j=(i+1):N + [I score_mat(i,j)] = cond_indep_chisquare(i,j,[],data,'LRT',alpha,node_sizes); + end +end +sc2=score_mat; + + +[tmp ordre]=sort(-score_mat(:)); +ordre2=ordre(find(-tmp>alpha)); +[II JJ]=ind2sub([N N],ordre2); + +pointer=1 ; +fini=length(II); + +% 3. +edges=2; +for pointer=1:min(2,fini), + %fprintf('%d-%d\n',II(pointer),JJ(pointer)); + G(II(pointer),JJ(pointer))=1; + G(JJ(pointer),II(pointer))=1; + score_mat(II(pointer),JJ(pointer))=-inf; +end + +pointer=min(2,fini); +arret=0; + +while pointer<fini & ~arret + % 4. + pointer=pointer+1; + C = ~reachability_graph(G); + if C(II(pointer),JJ(pointer)) + %fprintf('%d-%d\n',II(pointer),JJ(pointer)); + G(II(pointer),JJ(pointer))=1; + G(JJ(pointer),II(pointer))=1; + score_mat(II(pointer),JJ(pointer))=-inf; + edges=edges+1; + if edges==N-1 + arret=1; + end + end + + % 5. +end + +[tmp ordre]=sort(-score_mat(:)); +ordre2=ordre(find(-tmp>alpha)); +[II JJ]=ind2sub([N N],ordre2); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function G = phaseII(G1,data,node_sizes,alpha,II, JJ, score_mat) +% G = phaseII(G1,data,node_sizes,epsilon,II,JJ, score_mat) +% +% G1, II, JJ,score_mat are given by phaseI. +% data(i,m) is node i in case m. +% node_sizes is the vector of sizes ( default=max(data') ). +% alpha is the significant level for CI tests ( default=0.05 ). +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +st{1}='added'; +st{2}=''; +% 0. +[N m] = size(data); + +G=G1; + +% 6. +II=II(end:-1:1); +JJ=JJ(end:-1:1); + +pointer=length(II); + +while pointer>0 + % 7. + trysep = try_to_separate_A(G,II(pointer),JJ(pointer),data,alpha,node_sizes); + %fprintf('Try to separate %d and %d : %s\n',II(pointer),JJ(pointer),st{trysep+1}); + if ~trysep + G(II(pointer),JJ(pointer))=1; + G(JJ(pointer),II(pointer))=1; + %fprintf('%d-%d\n',II(pointer),JJ(pointer)); + end + % 8. + pointer=pointer-1; +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [G,U] = phaseIII(G1,data,node_sizes,alpha,s2,star) +% [G] = phaseIII(G1,data,node_sizes,alpha,s2,star) +% +% data(i,m) is node i in case m. +% if star~=0, use try_to_separate_B_star instead of try_to_separate_B (default star=1). +% G is the non-oriented graph and G1 the oriented result. +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +% 0. +if nargin < 2, disp('Not enough arguments');return; end +if nargin < 3, node_sizes=max(data'); end +if nargin < 4, alpha = 0.05; end +if nargin < 5, star=1; end +G=G1; +N=length(G); +% reachability_matrix of G +M = expm(full(G)) - eye(length(G)); M = (M>0); + +% 9. +fprintf('Thinning - separateA\n'); + +% Edges are examined in the inverse order of their Chi2 (or MI) score +s2(find(~G))=0; +[tmp ordre]=sort(s2(:)); +ordre2=ordre(find(tmp>0)); +[I J]=ind2sub([N N],ordre2); +ii=1:length(I); +for i=ii, + %fprintf('%d-%d : ',I(i),J(i)); + G(I(i),J(i))=0; + G(J(i),I(i))=0; + trysep = try_to_separate_A(G,I(i),J(i),data,alpha,node_sizes); + + if ~trysep, + G(I(i),J(i))=1; + G(J(i),I(i))=1; + %fprintf(' keep\n'); + %else + %fprintf('delete\n'); + end +end + +% 10. +fprintf('Thinning - separateB'); if star; fprintf('star'); end; fprintf('\n'); +s2(find(~G))=0; +[tmp ordre]=sort(s2(:)); +ordre2=ordre(find(tmp>0)); +[I J]=ind2sub([N N],ordre2); +ii=1:length(I); +for i=ii, + %fprintf('%d-%d : ',I(i),J(i)); + G(I(i),J(i))=0; G(J(i),I(i))=0; + if star==0 + trysep = try_to_separate_B(G,I(i),J(i),data,node_sizes,alpha); + else + trysep = try_to_separate_B_star(G,I(i),J(i),data,node_sizes,alpha); + end + if ~trysep, + G(I(i),J(i))=1; G(J(i),I(i))=1; + %fprintf(' keep\n'); + %else + %fprintf('delete\n'); + end +end + +%11. +fprintf('Thinning - orient_edges\n'); +U=G; +G=orient_edges(U,data,node_sizes,alpha); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [I, N1, N2] = try_to_separate_A(G,node1,node2,data,alpha,node_sizes) +% [I, N1, N2] = try_to_separate_A(G,node1,node2,data,alpha,node_sizess) +% +% G is the current partially directed graph. +% I is a boolean : I=1 <==> separated. +% N1 : neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2). +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +% 0. +if node1==node2 | length(G)<2 + disp('Error: Check your arguments in try_to_separate_A.');I=0; return +end +N=size(data,1); +if nargin==4, alpha=0.05; node_sizes=max(data'); end +if nargin==5, node_sizes=max(data'); end + +% 1. +N1=find(G(node1,:)==1);N01=N1; +GG1=G(setdiff(1:N,node1),setdiff(1:N,node1)); +node22=node2-(node2>node1); +% reachability_matrix of GG1 +M = expm(full(GG1)) - eye(length(GG1)); M = (M>0); +% N1 is the neighbors of node1 that are on the adjacency between node1 and node2 +for i=N1 + j=i-(i>node1); + if M(j,node22)~=1, N01=setdiff(N01,i); N1=setdiff(N1,i); end + % 2. + if ~G(node1,i), N1=setdiff(N1,i); end +end + +N2=find(G(node2,:)==1);N02=N2; +GG2=G(setdiff(1:N,node2),setdiff(1:N,node2)); +node12=node1-(node2<node1); +% reachability_matrix of GG2 +M = expm(full(GG2)) - eye(length(GG2)); M = (M>0); +% N2 is the neighbors of node2 that are on the adjacency between node1 and node2 +for i=N2 + j=i-(i>node2); + if M(node12,j)~=1, N02=setdiff(N02,i); N2=setdiff(N2,i); end + % 2. + if ~G(node2,i), N2=setdiff(N2,i); end +end + +%fprintf('%d : N1=',node1); fprintf('%d',N1); fprintf('\n'); +%fprintf('%d : N2=',node2); fprintf('%d',N2); fprintf('\n'); +% 3. +if length(N1)>length(N2), tmp=N1; N1=N2; N2=tmp; clear tmp, end +% 4. +C=N1; +for test=1:2 + if test==2, C=N2; end + % 5. + [I v1] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes); + if I, %fprintf('%d-%d separated (%2.5f) by C=',node1,node2,v1); fprintf('%d',C); fprintf('\n'); + return, + %else + %fprintf('%d-%d not separated (%2.5f) by C=',node1,node2,v1); fprintf('%d',C); fprintf('\n'); + end; + + + % 6. + step6=1; + while step6 + step6=0; + if length(C)>=1 + v=[]; + for i=C + Ci = setdiff(C,i); + [I(i) v(i)] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes); + end + [vm ind] = min(v); + + % 7. + if I(ind) + I = 1; return + else + if vm < v1 + v1 = vm; + C = setdiff(C,ind); + % goto step 6. + step6 = 1; + end + end + end + + % 8. + if test==2, I=0; return, end + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function I = try_to_separate_B(G,node1,node2,data,node_sizes,alpha,N1,N2) +% I = try_to_separate_B(G,node1,node2,data,node_sizes,epsilon,N1,N2) +% +% G is the current partially directed graph. +% data(i,m), node i in case m. +% node_sizes is the vector of size of the attributs in data ( default max(data') ) +% N1 (optionnal) is the neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2). +% I is a boolean : I+1 <==> separated. +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +% 0. +if node1==node2 | length(G)<2 + disp('Error: Verify your arguments in try_to_separate_B.');I=0; return +end +if nargin < 4, + disp('Error : not enougth arguments'); I=0; return; +end +N=size(data,1); + +% 1. +if nargin < 8 + N1=find(G(node1,:)==1); + GG1=G(setdiff(1:N,node1),setdiff(1:N,node1)); + node22=node2-(node2>node1); + M = expm(full(GG1)) - eye(length(GG1)); M = (M>0); + for i=N1 + j=i-(i>node1); + if M(j,node22)~=1, N1=setdiff(N1,i); end + end + N2=find(G(node2,:)==1); + GG2=G(setdiff(1:N,node2),setdiff(1:N,node2)); + node12=node1-(node2<node1); + M = expm(full(GG2)) - eye(length(GG2)); M = (M>0); + for i=N2 + j=i-(i>node2); + if M(node12,j)~=1, N2=setdiff(N2,i); end + end +end +if nargin < 6, alpha=0.05; end +if nargin < 5, node_sizes=max(data'); end +M = expm(full(G)) - eye(length(G)); M = (M>0); + +% 2. +N1b=[]; +for i=N1 + NN1=find(G(i,:)==1); + for j=NN1 + if M(i,j)~=1 & ~ismember(j,N1), N1b=union(N1b,j); end + end +end + +% 3. +N2b=[]; +for i=N2 + NN2=find(G(i,:)==1); + for j=NN2 + if M(i,j)~=1 & ~ismember(j,N2), N2b=union(N2b,j); end + end +end + +% 4. +if length(union(N1,N1b)) < length(union(N2,N2b)) + C=union(N1,N1b); +else + C=union(N2,N2b); +end + +% 5. +continu=1; +while continu + l=length(C); + %fprintf('%d',continu); + [I v] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes); + if I==1; return, elseif l<2, I=0; return, end + + % 6. + Cb=C; + for i=1:l + Ci=setdiff(C,C(i)); + [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes); + e = (v+1)/3; % e is a small value... + if I==1,return, elseif vi<v+e, Cb=setdiff(Cb,C(i)); end + end + + % 7. + if length(Cb) < l, C=Cb; else continu==0; I=0; return; end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function I = try_to_separate_B_star(G,node1,node2,data,node_sizes,alpha,N1,N2) +% I = try_to_separate_B_star(G,node1,node2,data,epsilon,node_sizess) +% +% G is the current partially directed graph. +% I is a boolean. +% N1 is the neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2). +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +% 0. +if node1==node2 | length(G)<2 + disp('Error: Verify your arguments in try_to_separate_B_star.');I=0; return +end +if nargin < 4, + disp('Error : not enougth arguments');I=0; return; +end +N=size(data,1); + +% 1. +if nargin < 8 + N1=find(G(node1,:)==1); + GG1=G(setdiff(1:N,node1),setdiff(1:N,node1)); + node22=node2-(node2>node1); + M = expm(full(GG1)) - eye(length(GG1)); M = (M>0); + for i=N1 + j=i-(i>node1); + if M(j,node22)~=1, N1=setdiff(N1,i); end + % 2. + if ~G(node1,i), N1=setdiff(N1,i); end + end + N2=find(G(node2,:)==1); + GG2=G(setdiff(1:N,node2),setdiff(1:N,node2)); + node12=node1-(node2<node1); + M = expm(full(GG2)) - eye(length(GG2)); M = (M>0); + for i=N2 + j=i-(i>node2); + if M(node12,j)~=1, N2=setdiff(N2,i); end + % 2. + if ~G(node2,i), N2=setdiff(N2,i); end + end +end +if nargin < 6, alpha=0.05; end +if nargin < 5, node_sizes=max(data'); end + +% 3. +if length(N1)>length(N2), tmp=N1; N1=N2; N2=tmp; end + +% 4. +C=N1; +l=length(C); +I=0; + +% 5. +for test=1:2 + continu=1; + %test + if test==2 & ~isempty(N2) + C=N2; l=length(C); IsInCi=zeros(1,l); IsInCi(l)=1; + else IsInCi=zeros(1,l); + end + s=ones(1,l); + % Pour tous les sous-ensemble Ci de C : + while continu & ~isempty(C) + Ci = setdiff(C.*IsInCi,0); + [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes); + if I, %fprintf('%d-%d separated (%2.5f) by C=',node1,node2,vi); fprintf('%d',Ci); fprintf('\n'); + return, + %else + %fprintf('%d-%d not separated (%2.5f) by C=',node1,node2,vi); fprintf('%d',Ci); fprintf('\n'); + end; + % if I, return, end + + if IsInCi==s, continu=0; + else + IsInCi(l)=IsInCi(l)+1; + + notOK=1; i=l; + while notOK & i>1 + if IsInCi(i)>s(i), IsInCi(i)=0; IsInCi(i-1)=IsInCi(i-1)+1; else notOK=0; end + i=i-1; + end + end + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function G1 = orient_edges(G,data,node_sizes,alpha) +% G1 = orient_edges(G,data,node_sizes) +% +% G is the partially directed graph. +% data(i,m), node i in case m. +% node_sizes is the vector of size of the attributs in data ( default max(data') ) +% +% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie" +% Jie Cheng, David Bell and Weird Liu. + +% 0. +if nargin < 4, alpha=0.05; end +if nargin < 3, node_sizes=max(data'); end +if nargin < 2, disp(' Require at least two arguments.'); return, end +N=length(G); +G1=G; + +% 1. +[Lnode1 Lnode2]=find(triu(1-triu(G),1)); %[Lnode1 Lnode2]=ind2sub([N N],find(triu(1-triu(G),1))); +for ii = 1:length(Lnode1), + node1=Lnode1(ii); + node2=Lnode2(ii); + %fprintf('%d %d\n',node1,node2); + N1 = find(G(node1,:)==1); + N2 = find(G(node2,:)==1); + if ~isempty(intersect(N1,N2)) + %fprintf('V1='); fprintf('%d ',N1); fprintf('\n'); + %fprintf('V2='); fprintf('%d ',N2); fprintf('\n'); + GG1 = G(setdiff(1:N,node1),setdiff(1:N,node1)); + node22 = node2-(node2>node1); + % reachability_matrix of GG1 + M = expm(full(GG1)) - eye(length(GG1)); M = (M>0); + % N1 is the neighbors of node1 that are on the adjacency between node1 and node2 + for i = N1 + j = i-(i>node1); + if M(j,node22)~=1, N1=setdiff(N1,i); + end + end + GG2 = G(setdiff(1:N,node2),setdiff(1:N,node2)); + node12 = node1-(node2<node1); + % reachability_matrix of GG2 + M = expm(full(GG2)) - eye(length(GG2)); M = (M>0); + % N2 is the neighbors of node2 that are on the adjacency between node1 and node2 + for i=N2 + j = i-(i>node2); + if M(node12,j)~=1, N2=setdiff(N2,i); end + end + %fprintf('%d %d\n',node1,node2); + %fprintf('N1='); fprintf('%d ',N1); fprintf('\n'); + %fprintf('N2='); fprintf('%d ',N2); fprintf('\n'); + + % 2. + M = expm(full(G)) - eye(length(G)); M = (M>0); + N1b=N1; + for i=N1 + NN1 = find(G(i,:)==1); + for j=NN1 + if M(i,j)~=1 & ~ismember(j,N1), N1b=union(N1b,j); end + end + end + %fprintf('N1''='); fprintf('%d ',N1b); fprintf('\n'); + + % 3. + N2b=N2; + for i=N2 + NN2 = find(G(i,:)==1); + for j=NN2 + if M(i,j)~=1 & ~ismember(j,N2), N2b=union(N2b,j); end + end + end + %fprintf('N2''='); fprintf('%d ',N2b); fprintf('\n'); + + % 4. + if length(N1b) < length(N2b) + C = N1b; + else + C = N2b; + end + %l=length(C); + %fprintf('C='); fprintf('%d ',C); fprintf('\n'); + + % 7. + step5=1; + while step5 + step5=0; + %fprintf('.'); + % 5. + l=length(C); + [I v] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes); + %fprintf('C='); fprintf('%d ',C); + %fprintf(': %d %2.5f\n',I,v); + + step8=0; + if I==1 & v~=0 % v < epsilon + step8=1; + else + if l==1 + G1(C,node1)=0; G1(C,node2)=0; + fprintf('%d -> %d <- %d\n',node1,C,node2); + step8=1; + end + end + %fprintf('%d\n',step8); + + % 6. + if ~step8 + Cb=C; + for i=1:l + Ci=setdiff(C,C(i)); + [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes); + % e = (v+1)/3; % e is a small value... + if I==1 % vi < v+e + Cb=setdiff(Cb,C(i)); + if ismember(C(i),N1) & ismember(C(i),N2) + G1(C(i),node1)=0; G1(C(i),node2)=0; + fprintf('%d -> %d <- %d\n',node1,C(i),node2); + end + if I==1 % vi < epsilon + step8=1; + end + end + end % for + end % if + + % 7. + if ~step8 + if length(Cb) < length(C), C=Cb; end + if length(C) > 0, step5==1; end + end + end % while step5 + % step8 : passer � la paire de noeud suivant + %fprintf('\n'); + %else + %fprintf(' No common neighbor\n'); + end % if +end % for + +% 11. +step9=0; +fprintf('Infering directions '); +test = pdag_to_dag(G1); +while ~isdag(test) %& step9<N + step9=step9+1; + %if ~isempty(test) + %fprintf('.'); + % 9. + for a=1:N, for b=1:N, for c=1:N, + if a~=b & b~=c & c~=a + if G1(a,b)==1 & G1(b,a)==0 + %fprintf('%d -> %d ... \n',a,b); + if G1(b,c)==1 & G1(c,b)==1 + if G1(a,c)+G1(c,a)==0, + G1(c,b)=0; + fprintf('%d -> %d (9)\n',b,c); + end + end + end + end + end, end, end + + % 10. + for a=1:N-1, for b=a+1:N + if G1(a,b)==1 & G1(b,a)==1 + GGG1=xor(G1,G1'); % matrice des arcs orient�s de G1 + M = expm(full(GGG1)) - eye(length(GGG1)); M = (M>0); + if M(a,b)==1, + G1(b,a)=0; + fprintf('%d -> %d (10)\n',a,b); + end + end + end, end + + test = pdag_to_dag(G1); + if isempty(test), G1=return_one_edge(G1); end + %else + % G1, return, + %end +end % while 11. +G1 = pdag_to_dag(G1); +fprintf('%d boucles\n',step9); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function 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 = b & find(sum(sum(eye(length(G)).*M))); % is there no cycle ? +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function G=return_one_edge(G) +N=length(G); fam=[]; node2=[]; +permnode = randperm(5); i=0; %fprintf('rev\n'); +while isempty(fam) | isempty(node2) | i==N + i=i+1; + node = permnode(i); %node = ceil(rand(1)*N); + fam=find(G(node,:)==1); + par=find(G(:,node)==1); + fam = myunion(fam, par); + if ~isempty(fam) + node2 = fam(ceil(rand(1)*length(fam))); + par2=find(G(:,node)==1); + if ~isempty(intersect(par2, node)), node2=[]; end + end +end +if isempty(myintersect(node, par2)), + G(node, node2)=0; + G(node2, node)=1; + fprintf('%d -> %d (Rev)\n',node, node2); +else + G(node, node2)=1; + G(node2, node)=0; + fprintf('%d -> %d (Rev)\n',node2, node); +end |