pdag=mgraph_meek4ruls(G,pdag)

function pdag=mgraph_meek4ruls(G,pdag)
%Input G: is undirected adjecent matrix, pdag is the pattern
%Output: pdag is complete directed graph
%
% rule 1 step D.1
  [A,B] = find(pdag==-1); % a -> b
 
  for i=1:length(A)
    a = A(i); b = B(i);
    C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a 
    %check this line oct 09.2003
    if (~isempty(C) & b<C)
      pdag(b,C) = -1; pdag(C,b) = 0;
      %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
    end
  end
  A=[]; B=[]; C=[];
  % rule 2 Step D.2 , Rule 2. 3 may unstable according to Tetrad4.2
  [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)' ) ;
      if (a<b & (pdag(a,b)==1 & pdag(b,a)==1))
        pdag(a,b) = -1; pdag(b,a) = 0;
      %  fprintf('rule 2: %d -> %d\n', a, b);
      end
    end
  end
  A=[]; B=[];
  % % rule 3
   [A,B] = find(pdag==1); % a-b
   for i=1:length(A)
     clear G2  
     a = A(i); b = B(i);
     C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
     
     % C contains nodes c s.t. a-c->ba
     G2 = setdiag(G(C, C), 1);
    
     if (any(G2(:)==0) & a<b) % there are 2 different non adjacent elements of C
       pdag(a,b) = -1; pdag(b,a) = 0;
      %fprintf('rule 3: %d -> %d\n', a, b);
     end
   end
 
 %rule 4 from Meek(95), this rule may unstable
 %added wang 10.2003
 %clear A B C1 C2 G2 
 %close all
 %figure
 %pdag=[0 1 0 1; 1 0 0 0 ; -1 -1 0 0 ; 1 0 -1 0]
 %draw_graph(abs(pdag));
 %G=[0 1 1 1; 1 0 1 0 ; 1 1 0 1 ; 1 0 1 0]
 
 [A, B]=find((G)==1);
  for i=1:length(A)
     C1=[]; C2=[]; C12=[]; C1_pa=[]; C1_ch=[]; G2=[];
     a=A(i); b=B(i);
     C1_pa =find( (pdag(:,a)==-1)' ); %a's parients directed neighbor
     C1_ch=find((pdag(a,:)==-1)); %a's child directed neoghbor
     
     C1=find( G(a,:)==1) ;%a's undirected neighbor 
     C2=find( G(b,:)==1);  %b's undirected neighbor
     C12=intersect(C1,C2); % a.b common neighbor
     G2 = setdiag(G(C12, C12), 1);
     %[~isempty(C12), any(G2(:)==0) ,~isempty(intersect(C1_ch,C12)), ~isempty(intersect(C1_pa,C12))]
     %pause
     if ( (~isempty(C12) & any(G2(:)==0) ) & ( ~isempty(intersect(C1_ch,C12)) & ~isempty(intersect(C1_pa,C12)) ) )
         %non adjacent element in C12 , orign a's child with b as same direction of a's parients
         need_origned_node=intersect(C1_ch,C2);
         if ~isempty(need_origned_node)
            for jj=1:length(need_origned_node)
                if (pdag(b,need_origned_node(jj))==1 & pdag(need_origned_node(jj),b)==1)
                %check once more
                    if b< need_origned_node
                        pdag(b,need_origned_node)=-1; pdag(need_origned_node,b) = 0;
                        %fprintf('rule 4: %d -> %d \n', b,need_origned_node(jj));
                    end
                   % a ,b
                   % need_origned_node
                   % warndlg('Error in Rule 4; PLEAE CHECK');
                   % pdag
                   % G
                   % break;
                end
            end %end for jj
        end     
     end
 end