edmonds.m

function[ED]=edmonds(V,E)%
%input is a directed graph, 
%V- Set of vertices [v1, v2,v3....]
%E- Set of edges is [ v1 v2 weight(v1,v2); ...] format

% Author: Ashish Choudhary, Ph.D
%         Pharmaceutical genomics division, TGEN.
%         13208,E Shea Blvd, Scottsdale, AZ 85259.
% Note/Disclaimer: This code is for academic purposes only.
% The implementation is derived from the book by Alan Gibbons.

%initialization
ED(1).BV=[];% Bucket of vertices
ED(1).BE=[];% Bucket of Edges

ED(1).V=V;
ED(1).E=E;

CURRENT_i=1;

while(1)    % breaking condition later
    
    CURRENT_i
    V=ED(CURRENT_i).V
    E=ED(CURRENT_i).E
     
    VERTICES_NOT_IN_BV=setdiff(V,ED(CURRENT_i).BV);
    
    if (numel(VERTICES_NOT_IN_BV)==0)
        break; %first phase
    end
    
    %let us add the first such
    VERTEX_ADDED=VERTICES_NOT_IN_BV(1);
    ED(CURRENT_i).BV=[ED(CURRENT_i).BV,VERTEX_ADDED];%
    
    %now check if largest incoming edge has a positive value
        
    [EDGE_VALUE,EDGE_ADDED]=index_of_max_value_incoming_edge(E,VERTEX_ADDED);
    
    if EDGE_VALUE>0 %dont do anything otherwise
        
    %upon adding check if there is a cicuit. 
    %If so, we will need to relabel everything  
    [dist,path]=iscycle([E(ED(CURRENT_i).BE,:)],E(EDGE_ADDED,1),E(EDGE_ADDED,2))
    
  
    ED(CURRENT_i).VERTICESINCKT=path;
    %now if the path was of finite length this 
    
     % now  adding to edge buckets
    ED(CURRENT_i).BE=[ED(CURRENT_i).BE,EDGE_ADDED];
    if dist<Inf              
        [GSTR,MAPVERT,MAPEDGE]=relabelgraph(ED(CURRENT_i));
        ED(CURRENT_i).MAPPINGVERT=MAPVERT;
        ED(CURRENT_i).MAPPINGEDGE=MAPEDGE;
        
        CURRENT_i=CURRENT_i+1
        ED(CURRENT_i).BV=GSTR.BV;
        ED(CURRENT_i).BE=GSTR.BE;
        ED(CURRENT_i).V=GSTR.V;
        ED(CURRENT_i).E=GSTR.E;
        
    end 
   
    
    end% end of EDGE_VALUE>0
end
    
%And now the reconstruction phase
   
%TREEMAX=reconstruct(ED);



%%
function[FLAGS]=exists_incoming_edge(G,NODEARRAY)
%finds if the nodes listed in the array have an incoming edge

for i=1:numel(NODEARRAY)
    FLAGS(i)=ismember(NODEARRAY(i),G(:,2));
end



%%
function[EDGE_VALUE,EDGE_INDEX]=index_of_max_value_incoming_edge(G,VERTEX)
%first find incoming edges
if numel(G)==0
    EDGE_VALUE=-1;
    EDGE_INDEX=NaN;
    return;
end
%
INDICES=find(G(:,2)==VERTEX)
VALUES=(G(INDICES,3));
[EDGE_VALUE,LOC]=max(VALUES);
EDGE_INDEX=INDICES(LOC);


%%
function[dist,path]=iscycle(G,S,D)
   if size(G,1)>0
       MAXN=max([S,D,max( unique(G(:,1:2)) )]);
       G=[G;MAXN,MAXN,0];
    DG = sparse(G(:,1),G(:,2),G(:,3));
    
        
  
    [dist,path]=graphshortestpath(DG,S,D);
    %if there is no path from S to D try D to S
    if dist==Inf
        [dist,path]=graphshortestpath(DG,D,S);
    end 
   else
     dist=Inf;
     path=[];
   end