function [maxFlow, Flow]= FordFulkersonShow(EdgeList,EdgeWgts, Flow)
%[maxFlow, Flow]= FordFulkersonShow(EdgeList,EdgeWgts,Flow)
%program for applying the Ford-Fulkerson algorithm for finding the maximal
%flow in a directed wieghted nework.  The program follows Algorithm 9.13 in
%the book; an efficient implementation of Edmonds and Karp.  

%Inputs:  EdgeList = 2-row matrix of directed edges.  Vertices are numbered from 1 (source)
%to n+2(sink) with the numbers in between being transit vertices.  
%EdgeWgts = column vector of corresponding weights (edge capacities).
%Flow = (optional input) initial flow; should be admissible.  If none is
%specified, the zero flow will be used initially.

%Outputs:  maxFlow = positive number representing maximum network flow.
%Flow = column vector for an admissible flow with total flow = maxFlow
%It is assumed that the edge weights are positive integers.
numVer = max(max(EdgeList)); %number of vertices in flow network
%Begin with zero flow if none is specified:
if nargin <3
    Flow =zeros(size(EdgeWgts));
end
%Assuming the inputted flow is admissible, we initialize maxFlow variable
%to be the value of the initial flow:
a = find(EdgeList(:,1)==1); % a = indices of edges flowing from s
maxFlow=sum(Flow(a)); 

%Vertex labelling scheme:
%The source vertex (vertex #1) will have label [0, +1, Inf]
%All other vertices will adhere to the following labelling scheme:
%label(i) = [j, +/-1, f],
%where j is the parent vertex from which vertex i was scanned,
%if the augmented graph edge that resulting in scanning vertex i from
%vertex j was a forward edge, the second component is +1, if it was a
%backward edge, the second component is -1.  The third component is the
%minimum of the third component of Label(j) and the excess capacity (if the
%scan used a foward edge), or the flow (if the scan used a backward edge)
%on the augmented digraph edge that led to Label(j).
Label(1,:) = [0 1 Inf];

%Initialize flag = 0; for main while loop.  Flag will be set = 1 when max flow is found. 
flag = 0;

%begin main while loop
iterCount = 1;
while flag == 0 %max flow not yet found
    %Step 2:  Initialize Vectors
    UnscannedVertices = [1];
    UnlabeledVertices = [2:numVer];
    %Step 3: Scan and Label
    while length(UnscannedVertices) > 0 & ismember(numVer,UnlabeledVertices)
        v = UnscannedVertices(1); %pop a vertex for scanning
        f = Label(v,3); %flow capacity up to v
        ULVertexIndex = 1;
        for w = UnlabeledVertices
            a = find(EdgeList(:,1)==v & EdgeList(:,2)==w); 
            %a = edge index of edge from v to (if this edge is in network,
            %otherwise a will be empty matrix
            b = find(EdgeList(:,1)==w & EdgeList(:,2)==v);
            if ~isempty(a) & EdgeWgts(a) > Flow(a) 
                Label(w,:) =[v  1  min(f,EdgeWgts(a) - Flow(a))];
                UnlabeledVertices(ULVertexIndex) = []; %delete w from unlabelled vertices
                ULVertexIndex = ULVertexIndex - 1;
                UnscannedVertices = [UnscannedVertices w]; %append w to unscanned vertices
            elseif  ~isempty(b) &  Flow(b) > 0 & v ~= 1 %& w ~= numVer  --don't need, loop will exit before t is scanned
                Label(w,:) =[v  -1  min(f,Flow(b))];
                UnlabeledVertices(ULVertexIndex) = []; %delete w from unlabelled vertices
                ULVertexIndex = ULVertexIndex - 1;
                UnscannedVertices = [UnscannedVertices w]; %append w to unscanned vertices
            end
            ULVertexIndex = ULVertexIndex + 1;
        end
        UnscannedVertices(1) = [];
    end
    if ismember(numVer,UnlabeledVertices) %sink vertex could not be labelled, therefore flow is maximum
        flag = 1;
    end
    %Step 4: Augmenting the flow (if flow is not max)
    if flag == 0
        %Flow augmentation
        vertex = numVer; %start with sink vertex and backtrack
        extraFlow = Label(vertex,3);
        while vertex > 1 %backtrack until we get to source
        if Label(vertex,2)>0 %forward edge leading to vertex
            %need edge number:
            %Forward edge
            a = find(EdgeList(:,1)==Label(vertex,1) & EdgeList(:,2)==vertex);
            Flow(a) = Flow(a) + extraFlow;
        else %backward edge going from vertex
            %need edge number:
            a = find(EdgeList(:,1)==vertex & EdgeList(:,2)==Label(vertex,1));
            Flow(a) = Flow(a) - extraFlow; 
        end
        vertex = Label(vertex,1);
        end
        maxFlow = maxFlow + extraFlow;
    end
    fprintf('After iteration #%d, we have:\r', iterCount)
    fprintf('Labels:\r')
    Label
    fprintf('Flow:\r')
    Flow
    fprintf('Value of flow:\r')
    maxFlow
    fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\r\r')
    clear Label
    Label(1,:) = [0 1 Inf];
    iterCount =iterCount + 1;
end