function [FLOW,Asort,C1,C2,VF,flg] = NetworkFlow(A)
%
% This code computes a max-flow for the given 
% Network.
% It is not optimized for efficiency.
%
% Written by Ming Gu for Math 170, April 2007
%
% On input:
%
% Each row of A has three numbers [i,j, cap]
% indicating an edge (i,j) with capacity cap>=0.
%
% Each edge (i,j) must be included in A twice, 
% one row for cap on edge (i,j), the other for 
% cap on (j,i), even when one of the caps might be 0.
%
% This code assumes that node 1 is the source
% and node n the sink, where n = max(A(:,2)).
%
%
% On output:
%
% FLOW is the max-flow, where 
% each row of FLOW has three numbers [i,j, flow]
% indicating an edge (i,j) with flow on the edge. 
%
% (C1, C2) are the min-cut on the Network.
%
% VF is the value of FLOW
%
% flg ~= 0 indicates code failure.
%
% Asort is a sorted version of the network.
%

%%
%% initialize
%%
flg = 0;
m = size(A,1);
n = max(A(:,2));
MaxFlow=0;
C1  = [];
C2  = [];
%%
%% sort the network.
%%
Asort = zeros(m,3);
Nw  = zeros(n,2);
alpha = 0;
for k=1:n
    link = find(A(:,1)==k);
    ln   = length(link); 
    Nw(k,:) = [alpha+1,alpha+ln];
    Asort(alpha+1:alpha+ln,:) = [k*ones(ln,1),A(link,2:3)];
    alpha   = alpha + ln;
end
MAX = max(Asort(:,3));
FLOW = [Asort(:,1:2), zeros(m,1)];
VF   = 0;
%
% step records the number of NetworkFlow steps 
% required to solve the problem.
%
step = 0;
%%
%% main loop
%%   
while (MaxFlow == 0)
  step = step + 1;
%
% find all unsaturated path reachable from source
%   
   C1 = [1];
   Path = [];
   C2 = [2:n]';
   unsat = 1;
   while (unsat == 1)
       unsat = 0;
%%
%% look for new unsaturated nodes in C2 that can be reached from C1
%%
       for k=1:length(C1)
           kc   = C1(k);
           [Link,INw,IC2]= intersect(Asort((Nw(kc,1):Nw(kc,2)),2),C2);
           if (length(Link) == 0)
               continue;
           end
           UnSat = find(Asort(Nw(kc,1)+INw-1,3)-FLOW(Nw(kc,1)+INw-1,3)>0 );
           if (length(UnSat) == 0)
               continue;
           end
           unsat = 1;
%%
%% nodes in UnSat are unsaturated
%%
           lUS   = length(UnSat);
           C1new = Asort(Nw(kc,1)+INw(UnSat)-1,2);  
           Path  = [Path; [ones(lUS,1)*kc, C1new(:)]];
           C1    = union(C1,C1new);
           C1    = C1(:);
           C2    = setdiff(C2,C1new);
           C2    = C2(:);
           break;
       end
       if (unsat == 0)
           MaxFlow = 1;
       end
%%
%% if sink is in C1, we can increase the flow
%%
       if (ismember(n,C1))
           MaxFlow = 0;
           unsat   = 0;
%
% construct an unsat path from source to sink.
%
           path    = [n];
           for kk=size(Path,1):-1:1
               if (Path(kk,2) == path(1))
                   path = [Path(kk,1);path];
               end
           end
           if (path(1)~=1)
              flg = -1;
              return;
           end
           INC = MAX*2;
           for kk=1:length(path)-1
                 kp = path(kk);
                 [Next,Ik,Ip] = intersect(Asort((Nw(kp,1):Nw(kp,2)),2),path(kk+1));
                 if (length(Next)==0)
                     flg = - 2;
                     return;
                 end
                 INC = min(INC,Asort(Nw(kp,1)+Ik-1,3)-FLOW(Nw(kp,1)+Ik-1,3));
           end
           if (INC <= 0)
               flg = - 3;
               return;
           end
%%
%% INC is the amount to increase in the flow. 
%% Adjust the network flow on the unsat path.
%%
           for kk=1:length(path)-1
                 kp = path(kk);
                 kp1= path(kk+1);
                 [Next,Ik,Ip] = intersect(Asort((Nw(kp,1):Nw(kp,2)),2),path(kk+1));
                 [Next1,Ik1,Ip1] = intersect(Asort((Nw(kp1,1):Nw(kp1,2)),2),path(kk));
                 FLOW(Nw(kp,1)+Ik-1,3)   = FLOW(Nw(kp,1)+Ik-1,3) + INC;
                 FLOW(Nw(kp1,1)+Ik1-1,3) = FLOW(Nw(kp1,1)+Ik1-1,3) - INC;
           end
           VF = VF + INC;
       end
   end
end
