%
% This is a sample program to solve the Transportation Problem:
%    min(sum_{i,j} C_{i,j}*X_{i,j})
%      s.t. sum_i X_{i,j} = d_j
%           sum_j X_{i,j} = s_i
%                 X_{i,j} >=0.
%
% where we assume C_{i,j}>=0, d_j>0, s_i > 0.
%
% Setup data. See NetworkFlow2 for more comments.
%
n = 6;
m = 10;
%
% generate random supply and demand, making sure sum(supply) = sum(demand).
%
S = rand(m,1);
D = rand(n,1);
D = D/sum(D)*sum(S);
M = 2*sum(S);
C = rand(m,n);
tol = 1e-12;
%
% initialize
%
U = min(C,[],2); U = U(:);
V = zeros(n,1);
step = 0;
Delta=0;
VF   = 0;
X    = zeros(m,n);
format short e;
%
% main loop
%
while (1)
%
% Display current results
%
   format long e
   disp([VF, Delta, sum(S.*U)+sum(D.*V)]);
   format short e;
%
% Setup Network
%
   LinkMat = [];
%
% use tol to avoid mistakes due to round-off.
%
   [I,J] = find(abs(repmat(U,1,n)+repmat(V',m,1)-C)<=tol);
   for k=1:m
      LinkMat = [LinkMat; [1 k+1 S(k)]; [k+1 1 0]];
   end
   for k=1:length(I)
       LinkMat = [LinkMat; [1+I(k),m+1+J(k),M];[m+1+J(k),1+I(k),0]];
   end
   for k=1:n
      LinkMat = [LinkMat; [k+m+1 m+n+2 D(k)]; [m+n+2,k+m+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);
      end
   end
   if (verify == 0)
       disp('NetworkFlow Failure');
       break;
   end
%
% Optimal solution is found if value of flow is sum(S).
%
   if (abs(VF-sum(S)) < tol)
      X   = zeros(m,n);
      INDX = find(FLOW(:,1)>1 & FLOW(:,1)< m+2 & ... 
                  FLOW(:,2)>m+1 & FLOW(:,1)< m+n+2 & ...
                  FLOW(:,3) > 0);
      for kk=1:length(INDX)
          X(FLOW(INDX(kk),1)-1, FLOW(INDX(kk),2)-(m+1)) = FLOW(INDX(kk),3);
      end
      break;
   end
%
% Update U and V to increase dual objective.
%
   I0 = find(C1>1 & C1 <= m + 1);
   I0 = C1(I0)-1;
   J0 = find(C1>m+1 & C1 <= m+n + 1);
   J0 = C1(J0)-m-1;
   J0c= setdiff((1:n)',J0);
   Delta = min(min(C(I0,J0c)-(repmat(U(I0),1,length(J0c))+repmat(V(J0c)',length(I0),1))));
   if (Delta <= 0)
      flg = -1;
      break;
   end
   U(I0) = U(I0)+Delta;
   V(J0) = V(J0)-Delta;
   step = step + 1;
end
