function [flow,labels] = LIS_maxflow(A, T, verbose)
% LIS_maxflow      Calls the max-flow (graph-cut) algorithm 
%
% Here the maxflowmex (maxflow) of Boykov-Kolmogorov's algorithm is preferred.
% implemented by Michael Rubinstein, WDI R&D and IDC, (2008)
% see: http://vision.csd.uwo.ca/code/
% and "An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy 
% Minimization in Vision.", by Yuri Boykov and Vladimir Kolmogorov.
% In IEEE Trans. on Pattern Analysis and Machine Intelligence (PAMI), Sep. 2004
%
% (I have unsucsessfully tried to use an alternative algorithm for max-flow
% min-cut, but not put much effort in it. It would, perhaps, be nice to have.
% This is the reason why most of this file is full of unused errornous code)
%
%  use: [flow,labels] = LIS_maxflow(A, T);
% ----------------------------------------------------------------------------
%  - A      the capacity network except source and sink, n-by-n sparse (double)
%           where n is number of nodes
%  - T      the capacities from source and to sink for each node, n-by-2 (double)
%           All nodes should be direct connected to source or sink, i.e. each
%           row of T should have (at least) one non-zero.
%  - flow   the value for max flow
%  - labels the labels (0 or 1) for each node
% ----------------------------------------------------------------------------
% ex: 
%  A = sparse([1 1 2 2], [3 4 3 4], [3 1 1 1], 4, 4);
%  T = [2,0; 3,0; 0,2; 0,3];
%  [flow,labels] = LIS_maxflow(A, T);

% another example
% A = [0,2,2,0,0,0; 2,0,2,0,0,0; 2,2,0,2,2,0; 0,0,2,0,2,2; 0,0,2,2,0,2; 0,0,0,2,2,0];
% A = sparse(A); T = [2,0;4,0;1,0;0,1;0,2;0,4];

mfile = 'LIS_maxflow';
if nargin < 3
    verbose = false;
end
flow = 0;

% path_maxflowmex = 'f:\matlab\maxflow-v3.01';    % expected location of mex-file
path_maxflowmex = my_matlab_path('maxflow');
if ~(exist('maxflowmex','file') == 3)
    if (exist(path_maxflowmex, 'dir') == 7)
        addpath(path_maxflowmex);
    else
        disp(['Catalog ',path_maxflowmex,' does not exist.']);
    end
end

if (exist('maxflowmex','file') == 3)
    % use the preferred implementation.
    %
    % Is this a correct description of labels ???
    % labels is 1 if flow from node to sink is smaller than capacity, else 0
    % The min-cut is then between source and node, node belong to sink
    % Also other nodes may belong to the 'sink'-group, what about them??
    [flow, labels] = maxflowmex(A, sparse(T));
    if verbose
        disp([mfile,': flow was ',int2str(flow)]);
    end

elseif exist('use_wrong_max_flow.txt','file')  
    % bruke heller denne: 
    % http://en.wikipedia.org/wiki/Relabel-to-front_algorithm
    % flow = 0;
    % labels = false(size(A,1),1);
    disp([mfile,': This probably makes a mess!']);
    n = size(A,1);   % number of nodes
    m = nnz(A);      % number of edges
    S = zeros(n,1);    % stored at each node
    maxS = 2*max( T(:,1)+reshape(full(sum(A)),n,1), T(:,2)+reshape(full(sum(A,2)),n,1) );
    maxloop = 200;
    [i,j,q] = find(A);
    E = [i,j,q,zeros(m,3)];   % one row for each edge
    flow_out_hist = zeros(maxloop,1);
    
    for loop = 1:maxloop
        F = zeros(size(T)); % flow
        % 1. fill the nodes
        F(:,1) = max(min(maxS-S, T(:,1)),0);
        S = S + F(:,1);
        % 2. level the nodes 
        E(:,6) = 0;            % the actual flow
        E(:,5) = S(E(:,1)) - S(E(:,2));       % E(:,5) is potential over edge
        E(:,4) = min(E(:,3),max(E(:,5),0));   % E(:,4) is maximal flow in edge
        %                                    restricted by potential and capasity
        for i = find(E(:,4)>0)'    % now, let it flow
            u = E(i,1); v = E(i,2);
            fi = min([S(u), maxS(v)-S(v), E(i,4), abs(S(v)-S(u))/2]); % 
            E(i,6) = fi;             % the actual flow
            S(u) = S(u) - fi;
            S(v) = S(v) + fi;
        end
        % 3. empty the nodes
        F(:,2) = min(S, T(:,2));
        S = S - F(:,2);
        % 4. report
        flow_out_hist(loop) = sum(F(:,2));
%         disp([' loop ',int2str(loop),', sum(S) = ',num2str(sum(S)), ...
%               ', ',int2str(nnz(E(:,3)==E(:,6))),' edges with max-flow', ...
%               ', flow in = ',num2str(sum(F(:,1))), ...
%               ', flow out = ',num2str(flow_out_hist(loop))]);
        if loop > 12
            if sum(flow_out_hist(loop-(5:9))) >= sum(flow_out_hist(loop-(0:4)))
                flow_out_hist = flow_out_hist(1:loop);
                break;
            end
        end
    end
    labels = false(n,1);
    for e = find(E(:,6)>0)';  % all edges with flow
        A(E(e,1),E(e,2)) = A(E(e,1),E(e,2)) - E(e,6);
        A(E(e,2),E(e,1)) = 0;
    end
    q = zeros(n,1); i1 = 1;
    u = find(T(:,2) > F(:,2));
    q(1:numel(u)) = u; 
    i2 = numel(u)+1;
    while (i1 < i2);
        u = q(i1);
        if ~labels(u)
            labels(u) = true;
            for v = (find(full(A(:,u))))'
                q(i2) = v;
                i2 = i2+1;
            end
        end
        i1 = i1+1;
    end
    
    % [flow, labels] = maxflowmex(A, sparse(T));
    flow = floor(flow_out_hist(loop)); 
    disp([mfile,': flow was ',int2str(flow)]);
    % figure(3); plot(flow_out_hist);
    
elseif exist('not_use_graphmaxflow','file')  
    % but this can not be used here because there often are reciprocal edges. 
    % use it anyway, but first direct all edges! (which is wrong!!!)
    disp([mfile,': try graphmaxflow (in Matlab bioinfo toolbox).']);
    [num_nodes, n] = size(A);
    if (num_nodes ~= n); 
        error([mfile,': Input matrix A is not square.']);
    end
    disp(' try to direct the graph to use it in graphmaxflow.');
    [u,v] = find(A);
    for e = 1:numel(u)
        % edge e has nodes u(e) and v(e)
        %    netto flow into u    -  netto flow into v
        d = (T(u(e),1)-T(u(e),2)) - (T(v(e),1)-T(v(e),2));
        % disp([' e=',int2str(e),' u=',int2str(u(e)),' v=',int2str(v(e)),' d=',num2str(d)]);
        if (d > 0) 
            A(v(e),u(e)) = 0;   % and A(u,v) is nonzero -> A(v,u) = 0
        elseif (d == 0) && (u(e) < v(e))
            A(v(e),u(e)) = 0;   
        end
    end
    if (nnz(A) < numel(u))
        disp(' The max flow will (probably) not be found but a smaller flow.');
        disp([' A matrix had ',int2str(numel(u)),' nonzeros, but here only ',...
             int2str(nnz(A)),' are kept.']);
    end
    cm = [0, sparse(T(:,1)'), 0; sparse(n,1), sparse(A), sparse(T(:,2)); sparse(1,n+2)];
    [flow, F] = graphmaxflow(cm, 1, n+2);
    disp([mfile,': flow was ',int2str(flow)]);
    % usually much smaller than flow from maxflowmex since direction of flow 
    % on each edge was locked before this maxflow-algorithm.
    % labels should be true for the nodes connected to sink
    labels = false(n,1);
    cm = ((cm-F) > 0); % rest capacity > 0
    q = zeros(n,1); i1 = 1; 
    u = find(cm(:,n+2));
    q(1:numel(u)) = u; 
    i2 = numel(u)+1;
    while (i1 < i2)
        u = q(i1);
        if ~labels(u-1)
            labels(u-1) = true;
            for v = (find(cm(:,u)))';
                if ~labels(v-1)
                    q(i2) = v;
                    i2 = i2 + 1;
                end
            end
        end
        i1 = i1 + 1;
    end
    
else
    disp([mfile,': does not find mex-file (or alternatives) for maxflow.']);
    
end



% ************ 
% tilsvarende c-kode i http://en.wikipedia.org/wiki/Relabel-to-front_algorithm
% function [flow, F] = relabel_to_front(C, s, t)

return
