function [maxFlow,flowMatrix,cuts] = graphmaxflow(G,S,D,NVPArgs)
%GRAPHMAXFLOW calculates the maximum flow in a directed graph.
% 
% [M,F,K] = GRAPHMAXFLOW(G,S,D) Calculates the maximum flow of the directed
% graph G from node S to node D. G is square adjacency matrix that
% represents a directed graph with non reciprocal edges. Nonzero entries in
% G determine the capacity of the edges. M is the maximum flow and F is a
% sparse matrix with all the flow values for every edge. F(i,j) is the flow
% from node i to j. K is a logical row vector indicating the nodes
% connected to S after calculating the minimum between S and D. 
% 
% GRAPHMAXFLOW(...,'METHOD',METHOD) selects the algorithm to use, options
% are: 
%    'Edmonds'    - Edmonds and Karp algorithm. Implementation is based on
%                   a variation called the "labeling algorithm". Time
%                   complexity is O(n*e^2). 
%   ['Goldberg']  - Goldberg algorithm uses the generic method known as
%                   "preflow-push". Time complexity is O(n^2*sqrt(e)).
% 
% GRAPHMAXFLOW(...,'CAPACITY',C) provides custom capacities for the edges.
% C is a column vector with one entry for every edge in G, traversed
% column-wise. 
%
% Example:
%   % Create a directed graph with 6 nodes
%   cm = sparse([1 1 2 2 3 3 4 5],[2 3 4 5 4 5 6 6],[2 3 3 1 1 1 2 3],6,6)
%   % Call the maximum flow algorithm between 1 and 6
%   [M,F,K] = graphmaxflow(cm,1,6)
%   % View graph with original flows
%   h = view(biograph(cm,[],'ShowWeights','on'))
%   % View graph with actual flows
%   view(biograph(F,[],'ShowWeights','on'))
%   % Show in the original graph one solution of the mincut problem
%   set(h.Nodes(K(1,:)),'Color',[1 0 0])
%
% See also: GRAPHALLSHORTESTPATHS, GRAPHCONNCOMP, GRAPHISDAG,
% GRAPHISOMORPHISM, GRAPHISSPANTREE, GRAPHMINSPANTREE, GRAPHPRED2PATH,
% GRAPHSHORTESTPATH, GRAPHTOPOORDER, GRAPHTRAVERSE.
%
% References: 
%  [1]	J. Edmonds and R.M. Karp "Theoretical improvements in the
%       algorithmic efficiency for network flow problems" Journal of the
%       ACM, 19:248-264, 1972. 
%  [2]	A.V. Goldberg "A New Max-Flow Algorithm" MIT Technical report
%       MIT/LCS/TM-291, 1985. 

%   Copyright 2006-2021 The MathWorks, Inc.
arguments
    G {bioinfo.internal.ValidationHelper.mustBeAdjacencyGraph}
    S(1,1) double {mustBePositive, mustBeInteger}
    D(1,1) double {mustBePositive, mustBeInteger}
    NVPArgs.Method(1,1) string {bioinfo.internal.ValidationHelper.mustBeMemberCaseInsensitive(NVPArgs.Method, ["Edmonds", "Goldberg"])} = "Goldberg"
    NVPArgs.Capacity(:,1) {mustBeNumeric, mustBeReal} = []
end

% Map algorithms to graph implementations
if strcmpi(NVPArgs.Method, "Goldberg")
    algorithm = 'pushrelabel';
else 
    % Edmonds
    algorithm = 'augmentpath';
end

conversionArgs = {};
if ~isempty(NVPArgs.Capacity)
    if numel(NVPArgs.Capacity) ~= nnz(G)
        % This should be called weights, but because the argument is
        % capacity, check the error message here rather than in
        % biograph2matlab
        error(message('bioinfo:graphmaxflow:InvalidCapacities'));
    end
    conversionArgs = {'Weights', NVPArgs.Capacity};
end

G_graph = bioinfo.internal.biograph2matlab(G, conversionArgs{:});

[maxFlow, graphFlow, cutSource] = maxflow(G_graph, S, D, algorithm);

flowMatrix = adjacency(graphFlow, 'weighted');

cuts = false(1, numnodes(G_graph));
cuts(cutSource) = true;

