function [key, p] = DMST(G_cost)

%Find a Directed Minimum Spanning Tree (DMST) of a directed simpel graph
%
%Input:
%   G_cost: a MxM matrix containing the weights of the edges of the graph
%Note: 1). nodes of the graph is assumed to be from 1 to M and node M will be the root of the DMST
%      2). G_cost(i,j) represents the wieght of edge directed from node j to node i. Remember the graph is directed,
%          G_cost(i,j) is generally unequal to G_cost(j,i)
%      3). In constructing G_cost, if there is no edge directed from k to m, then assign G_cost(m,k) a large number
%          so that it will never be selected into the final DMST
%
%Output:
%   p: a 1xM vector containing the parent of each node from 1 to M
%   key: a 1xM vector with jth element being the weight of the edge directed from p(j) to j
%
%KEY CAREABOUTS:
%1. If there is an edge directed from node i to j, then i is a parent of j.
%2. This is a pretty old algorithm, though not seen in books about algorithms. Following are several good references:
%[1]. E. Lawler, "Combinatorial optimization: networks and matroids", Sauders College Publishing, 1976
%There is one short section about this algorithm. Very clear explanation and a good example.
%[2]. Y.J. Chu and T.H. Liu, "On the shortest arborescence of a directed graph", Science Sinica, v14, 1965, pp. 1396-1400
%The best referetence for this algorithm. Very good explanation, example and mathenatical proof
%[3]. J. Edmonds, ``Optimum branchings'', J. Research of the National Bureau of Standards, 71B, 1967, pp.233-240. 
%[4]. F. Bock, ``An algorithm to construct a minimum spanning tree in a directed network'', Developments in
%        Operations Research, Gordon and Breach, NY, 1971, pp. 29-44. 
%[5]. R. E. Tarjan, ``Finding Optimum Branchings'', Networks, v.7, 1977, pp.25-35. 
%3. Node(or vetex) M is the root of the final DMST. Its parent will be 0 and its cost is 0
%
%SEE ALSO MRPF_DMST.m, MRPF_DMST_CSE and mst.m

M=size(G_cost,1);
InWhichPN = zeros(1,M);
for i=1:M
    InWhichPN(i) = i;
end
NET=[]; node=[];
for j=1:M-1
    node.pn=j; %name of a node or a pseudo node, starting from M+1,...
    node.ov=j;  %node.ov stores the original vertices within the pseudonode
    %node.iov stores all the orignal vertices directed into node.pn
    %Note: in our case, all ov's in a pseudonode have the same iov,i.e, {1:M}\{all ov's in pn}
    %to determin which pseudonoe an original vertex belongs to, refer to InWhichPN
    if j==1
        node.iov=[2:M];
    else
        node.iov=[[1:j-1] [j+1:M]];
    end
    for i=1:length(node.iov)
        node.w(i)=G_cost(j, node.iov(i));
    end
    NET{j}=node;
end
noCycle=0; 
pseudonode=[]; %list the pseudonodes generated when contracting the network
%pseudonoe:pn, ov, iov; shall be updated right after it is generated (as a new pseudonoe)
%this removes the need of hierarchical expansion of the pseudonode during the final construction of the DMST
while ~noCycle
    %step I: for each node choose the minimum-weighted edge directed into that node, forming T(V,E)
    %Note: need to guarantee that node M will never be a node of NET
    %Note: a node could be an original vertex or a pseudonode
    %T: currect subgraph(the DMST algorithm stops when there is no cycle in T
    %T.V: node; T.E: node incident to T.V; T.cost: weight of the selected edge
    %NET: current graph and records all the nodes(except M) and their incident edges and their weights
    T=[];
    for j=1:length(NET)
        node = NET{j}; 
        if length(node.ov)==1 %the node is not a pesudonode
            T(j).V = node.ov;
            [Y I] = min(node.w);
            for k=1:length(node.iov)
                if node.iov(k)==M  break; end
            end
            if node.w(k) == Y
                T(j).E = M;
            else    
                T(j).E = node.iov(I);
            end
            T(j).cost=Y;
        else %the node is a pseudonode
            [Y I] = min(node.w, [], 2); 
            [Y1 I1] = min(Y);
            T(j).cost=Y1;
            for k=1:length(Y)
                IsM=0;
                if Y(k)==Y1
                    for p=1:length(node.iov)
                        if node.iov(p)==M break; end
                    end
                    if node.w(k,p)==Y1
                        T(j).V = node.ov(k); T(j).E = M; IsM=1;
                    end
                end
                if IsM  break; end
            end
            if ~IsM  T(j).V = node.ov(I1); T(j).E = node.iov(I(I1));  end
        end
    end
    %step II: 1). cyle detection on T; 2). if there is a cycle in T, contract the network or stop otherwise
    %1). cycle detection: detect all the cycles if any
    examined = zeros(1, length(T));
    W_tmp=[]; W_tmp_cost=[]; W_tmp_node=[];
    CYCLE=[]; noAdjNode=1;
    %NOTE: CYCLE.node;CYCLE.e stores the original edges in the cycle;CYCLE.cost stores the costs of the edges in the cycle
    while sum(examined) ~= length(T) %while not all edges in T have been examined
        if noAdjNode
            W_tmp=[]; W_tmp_cost=[]; W_tmp_node=[]; 
            for i=1:length(examined)
                if examined(i)==0 break; end
            end
            current_node = InWhichPN(T(i).V); examined(i)=1;
            W_tmp_node = [W_tmp_node current_node];
            W_tmp = [W_tmp; [T(i).V T(i).E]]; %add the newly found original edge
            W_tmp_cost = [W_tmp_cost T(i).cost];
            current_node=InWhichPN(T(i).E);
        end
        for i=1:length(T)
            noAdjNode=1;
            if (examined(i)==0) & (InWhichPN(T(i).V) == current_node) %there is an unexamined edge directed into "next_node"
                W_tmp = [W_tmp; [T(i).V T(i).E]]; 
                W_tmp_cost = [W_tmp_cost T(i).cost];
                W_tmp_node = [W_tmp_node current_node];
                current_node = InWhichPN(T(i).E); examined(i)=1; 
                noAdjNode=0; %find an adjacent node
                %check if a cycle is formed, i.e., check if current_node is in W_tmp_node
                noCycle=1;
                for k=1:length(W_tmp_node)
                    if current_node==W_tmp_node(k) %there is a cycle
                        noCycle=0;
                        p=length(CYCLE);
                        CYCLE(p+1).e = W_tmp(k:length(W_tmp_node), :);
                        CYCLE(p+1).cost = W_tmp_cost(k:length(W_tmp_node));
                        CYCLE(p+1).node = W_tmp_node(k:length(W_tmp_node));
                    end
                    if (~noCycle)  W_tmp=[]; W_tmp_cost=[]; W_tmp_node=[]; break; end
                end
            end
            if ~noAdjNode break; end
        end
    end
    if length(CYCLE)
        noCycle=0;
    else
        noCycle=1;
    end
    %2). if there is cycle(s) in T, contract the network accordingly 
    NET_new=[]; 
    IsNodeInCycle=zeros(1,length(NET)); %indicate if a node of NET is in a cycle or not
    %if a node is in a cycle, it will be removed from the NET
    if ~noCycle 
        p=length(pseudonode);
        for i=1:length(CYCLE) %for each cycle
            %form a new pseudonode
            pseudonode(p+i).pn = M+p+i;
            pseudonode(p+i).oe = CYCLE(i).e;
            pseudonode(p+i).ov = [];
            for k=1:length(CYCLE(i).node)%for each node in the cycle
                if CYCLE(i).node(k) > M %if it is a pseudonode, add the ov and oe in that pn into thie new pn
                    pseudonode(p+i).ov = [pseudonode(p+i).ov pseudonode(CYCLE(i).node(k)-M).ov];
                    pseudonode(p+i).oe = [pseudonode(p+i).oe; pseudonode(CYCLE(i).node(k)-M).oe];
                else
                    pseudonode(p+i).ov = [pseudonode(p+i).ov CYCLE(i).node(k)];
                end
            end
            %contract the network: partI. for the newly generated pseudonode 
            node=[]; node.pn = M+p+i; node.ov = pseudonode(p+i).ov;
            iov=[];
            for k=1:M
                if ~IsInSet(k, node.ov)
                    iov=[iov k];
                end
            end
            node.iov = iov; %put all the original vertices that are not in node.ov into node.iov
            %form a new node and update the weights of edges directed into this new node
            max_cycle_cost = max(CYCLE(i).cost);
            w_update = zeros(1, length(CYCLE(i).node));
            for k=1:length(CYCLE(i).node)
                w_update(k) = max_cycle_cost - CYCLE(i).cost(k);
            end
            for k=1:length(node.ov)
                for j=1:length(node.iov)
                    pp = InWhichPN(node.ov(k)); 
                    for q=1:length(NET) %ov is in the q-th node of the NET
                        if NET{q}.pn == pp  break; end
                    end
                    for qq=1:length(NET{q}.iov)
                        if NET{q}.iov(qq)==node.iov(j)  break; end
                    end
                    for qqq=1:length(NET{q}.ov)
                        if NET{q}.ov(qqq)==node.ov(k)  break; end
                    end
                    for r=1:length(CYCLE(i).node)
                        if pp == CYCLE(i).node(r)  break; end
                    end
                    node.w(k,j) = NET{q}.w(qqq,qq) + w_update(r);
                    IsNodeInCycle(q) = 1; %q-th node of the NET is in a cycle
                end
            end  
            NET_new{i} = node;
            %update the InWhichPN vector
            for k=1:length(pseudonode(p+i).ov)
                InWhichPN(pseudonode(p+i).ov(k)) = M+p+i;
            end
        end
        t = length(NET_new);
        for i=1:length(NET)
            if ~IsNodeInCycle(i)
                t = t+1; NET_new{t} = NET{i};
            end
        end
        NET=NET_new;
    end
end
%final step: expanding the pseudonodes to form a DMST from T and calculate the cost
%Note: the root of the DMST is the virtual vertex M
T_edge = []; %store all the original edges in T
for i=1:length(T)
    T_edge = [T_edge; [T(i).V T(i).E]];
end
PNinT=[];
for i=1:length(T)
    if (InWhichPN(T(i).V) ~= T(i).V) &(~IsInSet(InWhichPN(T(i).V), PNinT))
        PNinT = [PNinT InWhichPN(T(i).V)];
    end
    if (InWhichPN(T(i).E) ~= T(i).E) &(~IsInSet(InWhichPN(T(i).E), PNinT))
        PNinT = [PNinT InWhichPN(T(i).E)];  
    end
end
for i=1:length(PNinT)
    T_edge = [T_edge; pseudonode(PNinT(i)-M).oe];
end
p=zeros(1,M); key=zeros(1,M); 
ov_visited=[M];
examined = zeros(1,size(T_edge,1));
while sum(examined) ~= length(examined)
    noAdjNode=1;
    if noAdjNode
        for i=1:size(T_edge,1)
            if (~examined(i)) & (IsInSet(T_edge(i,2), ov_visited))
                break;
            end
        end
        ov_visited = [ov_visited T_edge(i,1)];
        p(T_edge(i,1)) = T_edge(i,2); 
        key(T_edge(i,1)) = G_cost(T_edge(i,1), T_edge(i,2));
        current_node = T_edge(i,2);
        for k=1:size(T_edge,1)
            if (~examined(k)) & (T_edge(k,1)==T_edge(i,1))
                examined(k)=1;
            end  
        end
    end
    for i=1:size(T_edge,1)
        noAdjNode=1;
        if (~examined(i)) & (T_edge(i,2) == current_node)
            ov_visited = [ov_visited T_edge(i,1)];
            p(T_edge(i,1)) = T_edge(i,2); 
            key(T_edge(i,1)) = G_cost(T_edge(i,1), T_edge(i,2));
            current_node = T_edge(i,2);
            noAdjNode=0; 
            for k=1:size(T_edge,1)
                if (~examined(k)) & (T_edge(k,1)==T_edge(i,1))
                    examined(k)=1;
                end  
            end            
        end
        if ~noAdjNode break; end
    end
end