function [f,TC,XFlg,nf,out] = mcnf(IJCUL,SCUL,varargin)
%MCNF Minimum cost network flow problem.
% [f,TC,XFlg,nf,out] = mcnf(IJCUL,SCUL,opt)
% IJCUL = [i j c u l],  arc data
%  SCUL = [s nc nu nl], node data
%     i = n-element vector of arc tails nodes, where n = number of arcs
%     j = n-element vector of arc head nodes
%     c = n-element vector of arc costs
%     u = n-element vector of arc upper bounds
%       = [Inf], default
%     l = n-element vector of arc lower bounds
%       = [0], default
%     s = m-element vector of net node supply, where m = number of nodes
%            s(i) > 0 => node i is supply node        (outflow > inflow)
%            s(i) = 0 => node i is transshipment node (outflow = inflow)
%            s(i) < 0 => node i is demand node        (outflow < inflow)
%    nc = m-element vector of node costs
%    nu = m-element vector of upper bounds on total node flow
%       = [Inf], default
%    nl = m-element vector of lower bounds on total node flow
%       = [-Inf], default
%     f = n-element vector of arc flows
%    nf = m-element vector of node flows
% Options:
%   opt = options structure, with fields:
%       .LP       = linear programming solver
%                 = 1, only use LPRSM
%                 = 2, only use LINPROG, if available
%                 = 3, (default) use LPRSM for n <= LrgProb, else LINPROG
%       .LrgProb  = number of arcs in large problem instance
%                 = 1000, default
%       .TolInt   = tolerance for rounding to integer using ISINT
%                 = [1e-6], default
%       .LPargOut = return c,Alt,blt,Aeq,beq,l,u struct w/o running LP
%                 = 1, do
%                 = 0, do not (default)
%       .Lprsm    = option structure for LPRSM
%       .Linprog  = option structure for LINPROG, with MCNF default:
%                 .Display = 'none'
%       = 'Field1',value1,'Field2',value2, ..., where multiple input
%         arguments or single cell array can be used instead of the full
%         options structure to change only selected fields
%   opt = MCNF('defaults'), to get defaults values for option structure
%    TC = total cost
%  XFlg = exitflag from LP solver
%   out = output structure, with fields:
%       .LP       = LP solver used
%       .Lmb      = set of Lagrangian multipliers from LINPROG
%       .Outrsm   = output structure from LPRSM

% Copyright (c) 1994-2006 by Michael G. Kay
% Matlog Version 9 13-Jan-2006 (http://www.ie.ncsu.edu/kay/matlog)

% Input Error Checking ****************************************************
linprogexist = exist('linprog','file');
if nargin == 3 && isstruct(varargin{:})     % Use opt only to check fields
   opt = struct('LP',[],'LrgProb',[],'TolInt',[],'LPargOut',[],...
      'Linprog',[],'Lprsm',[]);
else                                       % Set defaults for opt structure
   opt = struct('LP',3,...
      'LrgProb',1000,...
      'TolInt',1e-6,...
      'LPargOut',0,...
      'Linprog',[],...
      'Lprsm',lprsm('defaults'));
   if exist('linprog','file') == 2
      opt.Linprog = optimset(linprog('defaults'),'Display','none');
   end
   if nargin == 1 && strcmpi('defaults',IJCUL), f = opt; return, end
   error(nargchk(1,Inf,nargin))
end
if nargin > 2
   opt = optstructchk(opt,varargin);
   if ischar(opt), error(opt), end
end

[n,cIJCUL] = size(IJCUL);
if cIJCUL < 3 || cIJCUL > 5, error('IJCUL must be 3-5 column matrix.'), end
[i,j,c,u,l] = mat2vec(IJCUL);
I = list2incid([i j],1);
if isempty(u), u = Inf*ones(n,1); end
if isempty(l), l = zeros(n,1); end

if nargin < 2 || isempty(SCUL)
   SCUL = zeros(max(max([i abs(j)],[],1)),1);
end
[m,cSCUL] = size(SCUL);
if m == 1, SCUL = SCUL'; m = cSCUL; cSCUL = 1; end
[s,nc,nu,nl] = mat2vec(SCUL);
if isempty(nc), nc = zeros(m,1); end
if isempty(nu), nu = Inf*ones(m,1); end
if isempty(nl), nl = zeros(m,1); end

if any(j < 0) || any(i == j)
   error('All arcs must be directed without self-loops.')
elseif any(l > u) || any(l == Inf)
   error('Elements of "l" can not be greater than "u" or Inf.')
elseif size(I,1) ~= m
   error('Length of "s" must equal maximum node index in "i,j".')
elseif any(~isfinite(s))
   error('Elements of "s" must be finite.')
elseif sum(s(s >= 0)) < sum(s(s <= 0))
   error('Total supply cannot be less than total demand.')
elseif any(nl > nu)
   error('Elements of "nl" can not be greater than ''nu'' or Inf.')
elseif length(opt.LP(:)) ~= 1 || all(opt.LP ~= [1 2 3])
   error('Invalid "opt.LP" specified.')
elseif opt.LP == 2 && ~exist('linprog','file')
   error('LP Solver LINPROG not found -- use LPRSM instead.')
elseif length(opt.LrgProb(:)) ~= 1 || opt.LrgProb < 0
   error('Invalid "opt.LrgProb" specified.')
elseif length(opt.LPargOut(:)) ~= 1 || all(opt.LPargOut ~= [0 1])
   error('Invalid "opt.LPargOut" specified.')
elseif ~isempty(opt.Linprog) && ~isstruct(opt.Linprog)
   error('"opt.Linprog" must be a structure from LINPROG("defaults").')
elseif ~isempty(opt.Lprsm) && ~isstruct(opt.Lprsm)
   error('"opt.Lprsm" must be a structure from LPRSM("defaults").')
end
% End (Input Error Checking) **********************************************

[f,TC,XFlg,nf] = deal([]);
out = struct('LP',[],'Lmb',[],'Outrsm',[]);

% Node costs
c = c + (nc'*(I==1))';  % Add node cost nc(i) to arc [i j] cost

% Node bounds
isnb = nu < Inf | nl ~= 0;

% Augment network by adding nodes for nonzero finite node bounds
if any(isnb)
   a = (I(isnb,:) == 1) + 0;
   I(isnb,:) = I(isnb,:) - a;
   I = [I zeros(m,sum(isnb))];
   I(isnb,n+1:end) = eye(sum(isnb));
   I = [I; a -eye(sum(isnb))];
   
   s2 = zeros(sum(isnb),1);
   s2(s(isnb) < 0) = s(isnb & s < 0);
   s(isnb & s < 0) = 0;
   s = [s; s2];
   
   c = [c; zeros(sum(isnb),1)];
   u = [u; nu(isnb)];
   l = [l; nl(isnb)];
end

% Add dummy demand node if excess supply
excess = sum(s(s >= 0)) + sum(s(s <= 0));
idxsup = find(s > 0);
idxdem = find(s < 0);
idxdummy = [];
if excess > 0
   s = [s; -excess];
   idxdummy = length(s);
   idxdem = [idxdem; idxdummy];
   [i,j] = incid2list(I);
   i = [i; idxsup];
   j = [j; idxdummy*ones(length(idxsup),1)];
   c = [c; zeros(length(idxsup),1)];
   u = [u; Inf*ones(length(idxsup),1)];
   l = [l; zeros(length(idxsup),1)];
   I = list2incid([i j],1);
end

if opt.LPargOut == 1
   f = struct('c',c,'Alt',I,'blt',s,'Aeq',[],'beq',[],'LB',l,'UB',u);
   return
end

% Solve LP
if ((opt.LP == 3 && n > opt.LrgProb) || opt.LP == 2) && ...
      exist('linprog','file') == 2
   out.LP = 2;
   if opt.LPargOut == 1
      f = struct('c',c,'Alt',I,'blt',s,'Aeq',[],'beq',[],'LB',l,'UB',u);
      return
   end
   
   [f,TC,XFlg,a_,out.Lmb] = linprog(c,I,s,[],[],l,u,[],opt.Linprog);
   
else
   out.LP = 1;
   if opt.LPargOut == 1
      f = struct('c',c,'Alt',[],'blt',[],'Aeq',I(1:end-1,:),...
         'beq',s(1:end-1),'LB',l,'UB',u);
      return
   end
   
   [f,TC,XFlg,out.Outrsm] = lprsm(c,[],[],I(1:end-1,:),s(1:end-1),...
      l,u,[],opt.Lprsm);
end

% Solution feasibility
if nargout < 3
   if XFlg < 0
      error('LP solver failed (XFlg < 0).')
   elseif XFlg == 0
      warning(...
         'LP solver reached max number iter w/o converging (LP XFlg = 0).')
   end
elseif XFlg < 0
   return
end

% Make integer all flows within tolerance
f(isint(f,opt.TolInt)) = round(f(isint(f,opt.TolInt))); 
TC = c'*f;

% Add cost of demand nodes
s = s(1:m);
if any(s<0)
   TC = TC + -nc(s<0)'*s(s<0);
end

% Calc. node flows
if nargout > 3
   nf = (I == 1)*f;
   nf = nf(1:m);
   if ~isempty(idxsup) && ~isempty(idxdummy)
      nf(idxsup) = nf(idxsup) - f(j == idxdummy);
   end
   nf(s<0) = -s(s<0);
end

f = f(1:n);