%
% This is a sample program to solve the Optimal Flow problem.
%
% Setup data. See comments in NetworkFlow2 for more comments.
%
%
n = 3;
M = n + 1;
A = rand(n,n);
%
% initialize
%
U = max(A,[],2); U = U(:);
V = zeros(n,1);
step = 0;
Delta=0;
VF   = 0;
format short e;
%
% main loop
%
while (1)
%
% Display current results
%
   format long e
   disp([VF, Delta, sum(U)+sum(V)]);
   format short e;
%
% Setup Network
%
   LinkMat = [];
   [I,J] = find(repmat(U,1,n)+repmat(V',n,1)-A==0);
   for k=1:n
      LinkMat = [LinkMat; [1 k+1 1]; [k+1 1 0]];
   end
   for k=1:length(I)
       LinkMat = [LinkMat; [1+I(k),n+1+J(k),M];[n+1+J(k),1+I(k),0]];
   end
   for k=1:n
      LinkMat = [LinkMat; [k+n+1 2*n+2 1]; [2*n+2,k+n+1 0]];
   end
%
% Solve the NetworkFlow problem
%
   [FLOW,Asort,C1,C2,VF,flg] = NetworkFlow2(LinkMat);
   if (flg ~= 0)
       disp('NetworkFlow failure');
       break;
   end
%
% making sure the NetworkFlow problem was successfully solved.
%
   fcc = 0;
   for k=1:length(FLOW)
      verify = 1;
      if (ismember(FLOW(k,1),C1) & ismember(FLOW(k,2),C2))
          if (abs(FLOW(k,3)-Asort(k,3)) > 0)
              verify = 0;
              K      = k;
              break;
          end
          fcc = fcc + Asort(k,3);
%%          disp([k, fcc, FLOW(k,1), FLOW(k,2)]);
      end
   end
%%   disp([step,length(C1),VF,flg,verify,fcc]);
   if (verify == 0)
       disp('NetworkFlow Failure');
       break;
   end
%
% Optimal solution is found if value of flow is n.
%
   if (VF == n)
      X   = zeros(n,1);
      INDX = find(FLOW(:,1)>1 & FLOW(:,1)< n+2 & ... 
                  FLOW(:,2)>n+1 & FLOW(:,1)< 2*n+2 & ...
                  FLOW(:,3) == 1);
      if (length(INDX)~=n)
         disp('Something wrong with INDX');
         break;
      end
      X(FLOW(INDX,1)-1) = FLOW(INDX,2)-(n+1);
      break;
   end
%
% Update U and V to reduce dual objective.
%
   I0 = find(C1>1 & C1 <= n + 1);
   I0 = C1(I0)-1;
   J0 = find(C1>n+1 & C1 <= 2*n + 1);
   J0 = C1(J0)-n-1;
   if (length(I0)<=length(J0))
      flg = -2;
      break;
   end
   J0c   = setdiff((1:n)',J0);
   Delta = min(min(repmat(U(I0),1,length(J0c))+repmat(V(J0c)',length(I0),1)-A(I0,J0c)));
   if (Delta <= 0)
      flg = -1;
      break;
   end
   if (length(Delta)>0)
      U(I0) = U(I0)-Delta;
      V(J0) = V(J0)+Delta;
   end
   step = step + 1;
end
