learn_struct_bnpc.m.svn-base

bayesian network structrue learning matlab program

	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