Code covered by the BSD License

Edmonds algorithm

Ashish Choudhary (view profile)

An implementation of Edmond's algorithm to obtain the maximum spanning weight tree from a graph.

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

%now check if largest incoming edge has a positive value

if EDGE_VALUE>0 %dont do anything otherwise

%upon adding check if there is a cicuit.
%If so, we will need to relabel everything

ED(CURRENT_i).VERTICESINCKT=path;
%now if the path was of finite length this

% now  adding to edge buckets
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

```