image thumbnail

Ancestral polytree

by

 

01 Feb 2013 (Updated )

The code here provides the algorithm of learning ancestral polytree.

LAP_PolytreeRecover( CIMatric,PolytreeSkeletonMatric )
function [ PolytreeSkeletonMatric,MAPS,CINumber,AmbiguousArc,Success ]  = LAP_PolytreeRecover( CIMatric,PolytreeSkeletonMatric )
%% this function is created for learning ancestral polytree.


Dim = size( PolytreeSkeletonMatric,2 );

MAPS = zeros( 1,Dim );
ObscureArc = zeros( Dim );
VisitedPair = zeros( Dim );
CINumber = 0;
Success = 1;
% step 1: G1 contains the inner nodes in the network
G1Nodes = zeros(1,Dim);
t = 0;

for p = 1:Dim
    if sum( PolytreeSkeletonMatric( p,: ) ) > 1 % isequal( PolytreeSkeletonMatric( p,: ),zeros( 1,Dim ) ) == 1
        t = t + 1;
        G1Nodes( t ) = p; 
    end
end
G1Nodes = G1Nodes(1:t);


%h3 = view(biograph( G1Matric ))
L = size( G1Nodes,2 );
UnvisitedG1Nodes = zeros( 1,L );
VisitedNodes = zeros( 1,Dim );

t = 1;
%2 searching part
while isempty( find( UnvisitedG1Nodes == 0,1 ) ) == 0
    if t == 5
        %h3 = view(biograph( G1Matric ));
    end
    G1Matric = PolytreeSkeletonMatric( G1Nodes,G1Nodes );
    MaxEdge = -1 ;MaxNode = -1; Oriented = 0; 
    for p = 1:L
        if  UnvisitedG1Nodes( p ) == 0
            LocalEdge = 0 ;% sum( G1Matric( p,:) );
            for q = 1:L
                if G1Matric( p,q )==1 && G1Matric( q,p )==1
                    LocalEdge = LocalEdge + 1;
                end
            end
            if LocalEdge > MaxEdge
                MaxEdge = LocalEdge;
                MaxNode =  p;
            end
        end
    end

    % if MaxEdge == -1,break;end   % all nodes are visited
    UnvisitedG1Nodes( MaxNode ) = 1; 
    MaxNode = G1Nodes( MaxNode );
    VisitedNodes( MaxNode ) = 1;
    
    LocalEdge = 0;
    for p = 1:Dim
        if MaxNode ~= p && PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 1
            LocalEdge = LocalEdge + 1;
        end
    end
      
    while 1             % case a.  p -- V* -- q ==> p -> V* <- q
          Run = 1; p = 1;
          while p <= ( Dim-1 ) && Run == 1
              if  PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 1 
                  q = p+1;
                  while q <= Dim && Run == 1
                       if  PolytreeSkeletonMatric( MaxNode,q ) == 1 && PolytreeSkeletonMatric( q,MaxNode ) == 1 && VisitedPair( p,q ) == 0   
                          % if CINumber < 10
                          % pa = p
                          % qa = q
                           
                          % end
                          VisitedPair( p,q ) = 1;
                          VisitedPair( q,p ) = 1;
                          CINumber = CINumber + 1; 
                          if CIMatric( p,q ) == 1 
                             %pa=p
                             %qa=q
                                 
                             Run = 0;
                             Oriented = 1; 
                             PolytreeSkeletonMatric( MaxNode,q )= 0 ;
                             PolytreeSkeletonMatric( MaxNode,p )= 0 ;                                                                                                                           
                          end                    
                       end
                       q = q+ 1;
                  end
              end
              p = p + 1;
          end
          if p >= Dim ,break;end
    end
 % Oriented1 = Oriented  
%if t > 0
%   h3 = view(biograph( PolytreeSkeletonMatric([1,2,5,6,7],[1,2,5,6,7]) ))
%   t = t -1;
%end

   while 1         % case b.   p--A-->q   => p-->V* <-->q 
          Run = 1; p = 1;
          while p <= Dim && Run == 1
              if  PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 1 
                  q = 1;
                  while q <= Dim && Run == 1
                       if PolytreeSkeletonMatric( q,MaxNode ) == 0 && PolytreeSkeletonMatric( MaxNode,q ) == 1 && VisitedPair( p,q ) == 0   
                           pb=p;
                           qb=q;
                          CINumber = CINumber + 1; 
                          if CIMatric( p,q ) == 1 
                            % pc=p
                            % qc=q
                             Run = 0;
                             Oriented = 1;                        
                             PolytreeSkeletonMatric( MaxNode,p )= 0 ;                                                                                                                           
                             PolytreeSkeletonMatric( MaxNode,q )= 0 ;
                             ObscureArc( q,MaxNode ) = 1;ObscureArc( MaxNode,q ) = 1; 
                          end                    
                       end
                       q = q+ 1;
                  end
              end
              p = p + 1;
          end
          if p > Dim ,break;end
   end
   
    while 1      % case c.  p--A<--q
          Run = 1; p = 1;
          while p <= Dim && Run == 1
              if  PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 1 
                  q =1;
                  while q <= Dim && Run == 1
                       if  PolytreeSkeletonMatric( q,MaxNode ) == 1 && PolytreeSkeletonMatric( MaxNode,q ) == 0 && VisitedPair( p,q ) == 0 
                            pc=p;
                            qc=q;
                          CINumber = CINumber + 1; 
                          if CIMatric( p,q ) == 1 
                             %pb=p
                             %qb=q
                             Run = 0;
                             Oriented = 1;                        
                             PolytreeSkeletonMatric( MaxNode,p )= 0 ;                                                                                                                           
                          end                    
                       end
                       q = q + 1;
                  end
              end
              p = p + 1;
          end
          if p > Dim ,break;end
    end
 
    while 1          % case d.   p-->A-->q
          Run = 1; p = 1;
          while p <= Dim && Run == 1
              if  PolytreeSkeletonMatric( MaxNode,p ) == 0 && PolytreeSkeletonMatric( p,MaxNode ) == 1  
                  q = 1;
                  while q <= Dim && Run == 1
                       if   q ~= p && PolytreeSkeletonMatric( q,MaxNode ) == 0 && PolytreeSkeletonMatric( MaxNode,q ) == 1 && VisitedPair( p,q ) == 0 
                             pd=p;
                             qd=q;
                          CINumber = CINumber + 1; 
                          if CIMatric( p,q ) == 1 
                             %pc=p
                             %qc=q
                             Run = 0;
                             Oriented = 1;                                                                                                                                                                                
                             PolytreeSkeletonMatric( MaxNode,q )= 0 ;
                             ObscureArc( q,MaxNode ) = 1;ObscureArc( MaxNode,q ) = 1;                               
                          end                    
                       end
                       q = q+ 1;
                  end
              end
              p = p + 1;
          end
         if p > Dim ,break;end
     end
    
    
  %  Oriented2 = Oriented
     TempPolytreeSkeletonMatric = PolytreeSkeletonMatric;Num = 0;
     for p = 1:Dim
         if PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 0
             Num = Num + 1;
         end
     end
     
     while 1           % case e.   p<--A-->q 
          Run = 1; p = 1;
          while p <= Dim-1 && Run == 1
              if  PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 0  
                  q = p+1;
                  while q <= Dim && Run == 1
                       if   q ~= p  && PolytreeSkeletonMatric( MaxNode,q ) == 1 && PolytreeSkeletonMatric( q,MaxNode ) == 0 && VisitedPair( p,q ) == 0 
                             pe=p;
                             qe=q;
                          CINumber = CINumber + 1; 
                          if CIMatric( p,q ) == 1 
                             %pc=p
                             %qc=q
                             Run = 0;
                             Oriented = 1;                        
                             PolytreeSkeletonMatric( MaxNode,p )= 0 ;                                                                                                                           
                             PolytreeSkeletonMatric( MaxNode,q )= 0 ;
                             ObscureArc( q,MaxNode ) = 1;ObscureArc( MaxNode,q ) = 1; 
                             ObscureArc( p,MaxNode ) = 1;ObscureArc( MaxNode,p ) = 1; 
                          end                    
                       end
                       q = q+ 1;
                  end
              end
              p = p + 1;
          end
         if p == Dim  ,break;end
     end
     if mod(Num,2) == 1
          p = 1;
          while p<= Dim
              if  PolytreeSkeletonMatric( MaxNode,p ) == 1 && PolytreeSkeletonMatric( p,MaxNode ) == 0
                   q = 1; 
                  while q <= Dim && Run == 1
                       if p~=q && TempPolytreeSkeletonMatric( MaxNode,q ) == 1 && TempPolytreeSkeletonMatric( q,MaxNode ) == 0 && VisitedPair( p,q ) == 0 
                         %  if CINumber < Control
                             pe=p;
                             qe=q;
                          % end
                          CINumber = CINumber + 1; 
                          if CIMatric( p,q ) == 1 
                             %pc=p
                             %qc=q
                             Run = 0;
                             Oriented = 1;
                             PolytreeSkeletonMatric( MaxNode,p )= 0 ;                                                                                                                           
                             ObscureArc( p,MaxNode ) = 1;ObscureArc( MaxNode,p ) = 1; 
                          end                    
                       end
                       q = q+ 1;
                  end
                  break;
              end
              p = p + 1;
          end
     end

% case e , if A->B-C, then orient A->B->C     
     if Oriented == 1
          MAPS( MaxNode ) = 1;
     end
    p = 1;
    while p <= Dim   
        if  PolytreeSkeletonMatric( p,MaxNode ) == 1 && PolytreeSkeletonMatric( MaxNode,p ) == 0     
            break;
        end
        p = p + 1;
    end
    if p<= Dim
       for q = 1:Dim
            if PolytreeSkeletonMatric( q,MaxNode )==1 && PolytreeSkeletonMatric( MaxNode,q ) == 1
                PolytreeSkeletonMatric( q,MaxNode )=0;
            end
        end
    end         

end          
   
   % Oriented3 = Oriented 
   % find( MAPS == 1) 
   for p = 1:Dim
       if MAPS( p ) == 1
           PolytreeSkeletonMatric = BFS_LocalPolytree( PolytreeSkeletonMatric,p );
       end
   end
  % find( MAPS == 1)
   for p = 1 : Dim
       if MAPS( p )==1
          for q = 1:Dim
              if MAPS(q)==1 && TestStrictPathInGraph_Set( PolytreeSkeletonMatric,p,q ) == 1
                  MAPS( q ) = 0;
              end           
          end
       end
   end
 % find( MAPS == 1) 
 % h3 = view(biograph( PolytreeSkeletonMatric ))
   
t = 1;AmbiguousArc = 0;
for p = 1 : Dim-1
    for q = (1+p) : Dim
        if ObscureArc( p,q )==1
           AmbiguousArc(t,1) = p; %#ok<AGROW>
           AmbiguousArc(t,2) = q; %#ok<AGROW>
           PolytreeSkeletonMatric( p,q ) = 1;
           PolytreeSkeletonMatric( q,p ) = 1;
           t = t + 1;
        end
    end
end

MAPS = find( MAPS==1 );
for p = 1:Dim-1
    for q = p+1:Dim
        if PolytreeSkeletonMatric( p,q ) == 1 && PolytreeSkeletonMatric( q,p ) == 1
            Success = 0;
            return
        end
    end
end
   
end

Contact us