function [x,xr,TC,nevals,XFlg] = milp(c,A,b,vlb,vub,N,R,varargin)
//MILP Mixed-integer linear programming.
//Branch Bound Algorithm.
// [y,x,FVAL,nevals,XFlg] = milp(c,A,b,vlb,vub,N,R,opt)
// //        
// //            min c'*X    subject to:   A*X <= | = b 
// //            x real
// //            y integer
// //            X=[y;x]
// //   Last R variables not integer.
// //   VLB & VUB  defines a set of lower and upper
// //   bounds on the design variables, X, so that the solution is in
// //   the range VLB <= X <= VUB.  Use empty matrices for VLB and VUB
// //   if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; 
// //   set VUB(i) = Inf if X(i) is unbounded above.
//       N = First N constraints equality constraints (=0, default)
//       R = Last R variables not integer constrained (=0, default)
//    varargin define a options vector.
// Options:
//     opt = options structure, with fields:
//         .Display  =-1, no warning messages
//                   = 0, initial LP message (default)
//                   = 1, show steps of B&B evalutation
//                   = 2, also show TolFun changes
//         .Martin   = 0, don't use Martin cuts
//                   = 1, use Martin cuts if possible (default)
//         .SOS      = Special Ordered Sets of type one (SOS1)
//                   = 0, don't use (default)
//                   = 1, use if possible (and don't use Martin cuts)
//         .RCthresh = reduced cost threshold percent of obj for Martin cuts
//                   = 0.05 (default)
//         .TolInt   = integer tolerance for ISINT
//                   = [1e-8], default
//         .ExTolFun = opt.LPopt.TolFun expansion factor for LINPROG feasible
//                   = 10 (default)
//         .LPopt    = 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 = MILP('defaults'), to get defaults values for option structure
//
// OUTPUTS:
//  y returns  the value of the integers variables
//  x returns  the value of the real variables
//  FVAL returns the value of the objective function
//  nevals = number of LP evaluations
//    XFlg = 2, initial LP solution feasible
//         = 1, solution found
//         =-1, no feasible integer solution found
//         =-2, initial LINPROG exitflag <= 0

// Copyright (c) 1994-2003 by Michael G. Kay
// Modified by Marcelo Escobar at Sep-2004
// Matlog Version 7 02-Sep-2003

// Input Error Checking ******************************************************
if nargin == 8 && isstruct(varargin{:})     // Use opt only to check fields
   opt = struct('Display',[],'Martin',[],'SOS',[],'RCthresh',[],...
      'TolInt',[],'ExTolFun',[],'LPopt',[]);
else                                       // Set defaults for opt structure
   opt = struct('Display',0,...
      'Martin',1,...
      'SOS',0,...
      'RCthresh',0.05,...
      'TolInt',0.01*sqrt(eps),...
      'ExTolFun',10,...
      'LPopt',optimset(optimset('linprog'),'Display','none'));
   if nargin == 1 && strcmpi('defaults',c), x = opt; return, end
   error(nargchk(5,Inf,nargin))
end
if nargin > 7
   opt = optstructchk(opt,varargin);
   if ischar(opt), error(opt), end
end

if nargin < 6 || isempty(N), N = 0; end
if nargin < 7 || isempty(R), R = 0; end

c = c(:)'; n = length(c);
b = b(:); m = length(b);
[mA,nA] = size(A);
vlb = vlb(:)';vub = vub(:)';

if mA ~= m || nA ~= n
   error('Dimensions of c and b not compatible with A.');
elseif length(vlb) ~= n || length(vub) ~= n
   error('Variable bounds not the same dimension as c.');
elseif length(N(:)) ~= 1 || N < 0 || N > m
   error('Number of equality constraints N not feasible.');
elseif length(R(:)) ~= 1 || R < 0 || R > n
   error('Number of non-integer-constrained variables R not feasible.');
elseif length(opt.Display(:)) ~= 1 || ~or(opt.Display == [-1 0 1 2])
   error('Invalid ''opt.Display'' specified.')
elseif length(opt.Martin(:)) ~= 1 || ~or(opt.Martin == [0 1])
   error('Invalid ''opt.Martin'' specified.')
elseif length(opt.SOS(:)) ~= 1 || ~or(opt.SOS == [0 1])
   error('Invalid ''opt.SOS'' specified.')
elseif length(opt.ExTolFun(:)) ~= 1
   error('Invalid ''opt.ExTolFun'' specified.')
end
// End (Input Error Checking) ************************************************

// Convert options to use Martin and SOS1
Ar=A;Br=b;
cr=c;lbr=vlb;ubr=vub;

xr=[];
useMartin = 0; useSOS1 = 0;
if opt.Martin == 1, useMartin = 1; end
if opt.SOS == 1, useSOS1 = 1; end
if useSOS1 == 1 & useMartin == 1, useMartin = 0; end  // No Martin if SOS

// Initial LP solution
[x,fval,how,out,lmb] = linprog(...
   c,A(N+1:end,:),b(N+1:end),A(1:N,:),b(1:N),vlb,vub,[],opt.LPopt);
nevals = 1;

// Case 1: No integer variables
if n == R
   TC = c*x;
   XFlg = 2;
   return;
end

// Case 2: Integer variables
I = n - R;
// x(isint(x,opt.TolInt)) = round(x(isint(x,opt.TolInt)));
y = x(1:I);//**** x=

// Case 2a: LP solution not OK
if how <= 0
   TC = c*x;
   XFlg = -2;
   return;
// Case 2b: Initial LP solution feasible
elseif all(isint(x,opt.TolInt))
   TC = c*x;
   if opt.Display > 0
      disp('Initial LP solution feasible.');
   end
   XFlg = 2;
   return;
end

// Case 2c: Initial LP solution not feasible
TC = Inf;
TC0 = c*x;
x0 = x;
chow = how;

vlbR = vlb(I+1:n);vlbI = vlb(1:I);
vubR = vub(I+1:n);vubI = vub(1:I);

// Check if SOS1 possible for binary variables
if useSOS1 && any(vlbI == 0 && vubI == 1) && or(b(1:N) == 1)
   varB = vlbI == 0 && vubI == 1;
   cstr = 1:N;cstr = cstr(b(cstr) == 1);
   cstr = cstr(mtlb_all(A(cstr,1:I)' == 0 || A(cstr,1:I)' == 1));
   if R > 1
      cstr = cstr(mtlb_all(A(cstr,I+1:n)' == 0));
   elseif R == 1
      cstr = cstr(A(cstr,I+1:n)' == 0);
   end
   idSOS1 = 0;
   SOS1 = zeros(1,I);
   for i = 1:length(cstr)
      nSOS1 = sum(A(cstr(i),1:I) && varB);
      if nSOS1 > 2			// SOS1 <= 2 same as reg. B&B
	 idSOS1 = idSOS1 + 1;
	 SOS1(A(cstr(i),1:I) && varB) = idSOS1*ones(1,nSOS1);
      end
   end
   if mtlb_all(SOS1 == 0)
      useSOS1 = 0;
   end
else
   useSOS1 = 0;
end

if opt.Display > 0
   fprintf('\nNODE PRNT   OBJ       BND       VAR\n');
   fprintf('//4d      //10.4f\n',0,TC0);
end

// Fix variables via Martin cuts
fix = 0;				// fix = 1 => 1st Martin cut
release = 0;				// release = 1 => 2nd Martin cut 
ixnofix = 1:n;
if useMartin
   // [x(1:I) lmb(m+1:m+I) lmb(m+I+1:m+2*I)]
   RCthresh = opt.RCthresh * abs(TC0);		// Reduced cost > RCthresh // of obj 
   ixfix = find(x == 0 && lmb.lower > RCthresh);//********** find(x(1:I
   nfix = length(ixfix);
   if nfix
      ixnofix(ixfix) = [Inf];
      fix = 1;
      TC_init = TC0;
      c_init = c;c(ixfix) = [Inf];
      A_init = A;A(:,ixfix) = [Inf];
      vlbI_init = vlbI;vlbI(ixfix) = [Inf];
      vubI_init = vubI;vubI(ixfix) = [Inf];
      I = I - nfix;
      x(ixfix) = [0];x0(ixfix) = [Inf];
      if opt.Display > 0
	  fprintf('Fixing');fprintf(' //d',ixfix);
	  fprintf(' with RC > //.4f\n',RCthresh);
      end
      if useSOS1
	 SOS1_init = SOS1;SOS1(ixfix) = [Inf];
	 if ~or(SOS1)
	    useSOS1 = -1;
	 end
      end
   end
end

// Initialize B&B tree "T"
T = [TC0 vlbI vubI 0 NaN 0];		// ... 0 NaN 0] => ... prnt ixx bndx]
if useSOS1 > 0
   T = [T SOS1];
end

// Main B&B tree loop
while ~isempty(T) || fix

   // Release fixed variables (2nd Martin cut) after solved for 1st Martin cut
   if isempty(T) && fix
      fix = 0;release = 1;
      T = [TC_init vlbI_init vubI_init 0 NaN 0];
      if useSOS1 ~= 0;
	 T = [T SOS1_init];
	 useSOS1 = 1;
      end
      Arelease = zeros(1,n);Arelease(ixfix) = -ones(1,nfix);
      A = [A_init;Arelease];b = [b;-1];
      c = c_init;I = I + nfix;
      xynofix = x;x = zeros(n,1);
      if sum(ixnofix)>=0
//           x(ixnofix) = xynofix;
      end //////// condicional
      ixnofix = 1:n;
      if opt.Display > 0
	 fprintf('Releasing fixed variables\n');
      end
   end
   
   // Node selection
   [tmp,ixT] = min(T(:,1));		// Selection a.k.a. "jumptracking"
   
   // Get selected node info and then remove from B&B tree
   cvlbI = T(ixT,2:I+1);
   cvubI = T(ixT,I+2:2*I+1);
   prnt = T(ixT,2*I+2);			// Parent of current node
   ixx = T(ixT,2*I+3);		        // Index of branching var, where 
                                        //    ixx == NaN => initial root node
					//               or releasing fix var
					//    ixx < 0 => LB, ixx > 0 => UB
					//    ixx == 0 => using SOS1 var
   bndx = T(ixT,2*I+4);			// LB or UB of ixx
   if useSOS1 > 0
      SOS1 = T(ixT,2*I+5:3*I+4);
      if ixx == 0
         ixSOS1 = find(SOS1 == -1);	// For display, current SOS1 var = -1
         SOS1(ixSOS1) = zeros(1,length(ixSOS1));
      end
   end
   T(ixT,:) = [];
   
   // Solve LP relaxation
 
   if ~isnan(ixx) || release		// Not initial root node
      [x0,fval,chow] = linprog(...
         c,A(N+1:end,:),b(N+1:end),A(1:N,:),b(1:N),...
         [cvlbI vlbR],[cvubI vubR],[],opt.LPopt);      // x = X0
      
      if chow <= 0 // Overly constrained: Loosen TolFun to find feas sol
         origTolFun = opt.LPopt.TolFun;
         opt.LPopt.TolFun = opt.ExTolFun * origTolFun;
         if opt.Display > 1
            fprintf('//36s Increasing TolFun for next node\n',' ');
         end
         [x0,fval,chow] = linprog(...
            c,A(N+1:end,:),b(N+1:end),A(1:N,:),b(1:N),...
            [cvlbI vlbR],[cvubI vubR],[],opt.LPopt);   // x0 = X0
         opt.LPopt.TolFun = origTolFun;
      end
      
      nevals = nevals + 1;
   end
   
   // Current solution
   
   x0(isint(x0,opt.TolInt)) = round(x0(isint(x0,opt.TolInt)));
   x0 = x0(1:I);
   TC0 = c(1:I)*x0;//****c puro
   
   if chow > 0
      // Current solution cannot be fathomed
      if TC0 < TC && ~mtlb_all(isint(x0,opt.TolInt))
	 if opt.Display > 0 && (~isnan(ixx) || release)
	    fprintf('//4d //4d //10.4f ',nevals-1,prnt,TC0);
	    if ixx < 0
	       fprintf('//4d LB of //3d\n',bndx,ixnofix(-ixx));
	    elseif ixx > 0
	       fprintf('//4d UB of //3d\n',bndx,ixnofix(ixx));
	    elseif ixx == 0
	       fprintf('//4d UB of',bndx);fprintf(' //3d',ixnofix(ixSOS1));
	       fprintf('\n');
	    end
	 end	 
	 
	 // Select branching variable
         ixfrc = find(~isint(x0,opt.TolInt));
	 [tmp,ixmxfrc] = max(min(ceil(x0(ixfrc)) - x0(ixfrc),...
	                         x0(ixfrc) - floor(x0(ixfrc))));
         ixx = ixfrc(ixmxfrc);		// Maximum integer infeasibility
         dvubI = cvubI;dvubI(ixx) = floor(x0(ixx));
         uvlbI = cvlbI;uvlbI(ixx) = ceil(x0(ixx));
	 
	 // Add to B&B tree
     
    if ~useSOS1			// No SOS1
       T = [T;TC0 cvlbI dvubI nevals-1 ixx dvubI(ixx);...
             TC0 uvlbI cvubI nevals-1 -ixx uvlbI(ixx)];
       else				// Use SOS1
          lSOS1 = SOS1;rSOS1 = SOS1;
          if SOS1(ixx)		// In SOS1
             ixSOS1 = find(SOS1 == SOS1(ixx)); // Find all SOS1 for ixx
             if length(ixSOS1) <= 2	// Use reg. bnds. when <= 2 in SOS1
                lSOS1(ixSOS1) = zeros(1,length(ixSOS1));rSOS1 = lSOS1;
             else			// Use left & right SOS1 bnds. with a
                // balanced no. of fractional vars
                ixfrc = ixSOS1(~isint(x0(ixSOS1)));nfrc = length(ixfrc);
                ixint = ixSOS1(isint(x0(ixSOS1))); nint = length(ixint);
                ixl = [ixfrc(1:floor(nfrc/2)) ixint(1:ceil(nint/2))];
                ixr = [ixfrc(floor(nfrc/2)+1:nfrc) ...
                      ixint(ceil(nint/2)+1:nint)];
                lSOS1(ixl) = -ones(1,length(ixl)); // Set to -1 for display
                rSOS1(ixr) = -ones(1,length(ixr));
                lvubI = cvubI;lvubI(ixl) = zeros(1,length(ixl));
                rvubI = cvubI;rvubI(ixr) = zeros(1,length(ixr));
                ixx = 0;		// Set to 0 to indicate using SOS1
             end
          end	    
          if ixx ~= 0			// SOS1 not used, but pass l&r SOS1
             T = [T;TC0 cvlbI dvubI nevals-1 ixx dvubI(ixx) lSOS1;...
                   TC0 uvlbI cvubI nevals-1 -ixx uvlbI(ixx) rSOS1];
             else			// Use ixx = 0 and bndx = 0
                T = [T;TC0 cvlbI lvubI nevals-1 0 0 lSOS1;...
                      TC0 cvlbI rvubI nevals-1 0 0 rSOS1];
                end
             end
             // Current solution new incumbent
      elseif TC0 < TC && mtlb_all(isint(x0,opt.TolInt));
         x = x0;
         TC = TC0;
         if opt.Display > 0
            fprintf('//4d //4d //10.4f ',nevals-1,prnt,TC);
            if ixx < 0
               fprintf('//4d LB of //3d',bndx,ixnofix(-ixx));
            elseif ixx > 0
               fprintf('//4d UB of //3d',bndx,ixnofix(ixx));
            elseif ixx == 0
               fprintf('//4d UB of',bndx);fprintf(' //3d',ixnofix(ixSOS1));
            end
            fprintf(' Incumbent\n');
         end	 
         if ~isempty(T)
            T = T(T(:,1) < TC,:);
         end
         // Current solution fathomed due to value dominance
      else
         if opt.Display > 0
            fprintf('//4d //4d //10.4f ',nevals-1,prnt,TC0);
            if ixx < 0
               fprintf('//4d LB of //3d',bndx,ixnofix(-ixx));
            elseif ixx > 0
               fprintf('//4d UB of //3d',bndx,ixnofix(ixx));
            elseif ixx == 0
               fprintf('//4d UB of',bndx);fprintf(' //3d',ixnofix(ixSOS1));
            end
            fprintf(' Fathomed: value dominance\n');
         end	 
      end
   // Current solution fathomed due to not OK LP relaxation
   else
      if opt.Display > 0
         fprintf('//4d //4d //10.4f ',nevals-1,prnt,TC0);
         if ixx < 0
            fprintf('//4d LB of //3d',bndx,ixnofix(-ixx));
         elseif ixx > 0
            fprintf('//4d UB of //3d',bndx,ixnofix(ixx));
         elseif ixx == 0
            fprintf('//4d UB of',bndx);fprintf(' //3d',ixnofix(ixSOS1));
         end
         fprintf(' Fathomed: infeasible\n');  // (' Fathomed: //d\n',chow)
      end 
   end
end // End B&B tree evaluation

// No incumbent found during B&B tree evaluation
if TC == Inf
   if opt.Display > -1
      disp('Warning: There is no feasible integer solution.')
   end
   TC = c*x;
   XFlg = -1;
else
   XFlg = 1; 
   
end

if R
Cr=ones(1,1:R);
N=N+1;
Ar=[cr;Ar];
Br=[TC;Br];
Arm=Ar(:,end-R+1);
Brm=Br-Ar(:,1:end-R)*x;
vlbr=lbr(end-R+1);
vubr=ubr(end-R+1);
[xr] = linprog(Cr,Arm(N+1:end,:),Brm(N+1:end),Arm(1:N,:),Brm(1:N),vlbr,vubr,[]);
end

endfunction

//SUBROTINAS:
//-------------------------------------------------------------------------
function y = isint(x,TolInt)
//ISINT True for integer elements (within tolerance).
//      y = isint(x,TolInt)
//        = abs(x-round(x)) < TolInt
// TolInt = integer tolerance
//        = [0.01*sqrt(eps)], default

// Copyright (c) 1994-2003 by Michael G. Kay
// Matlog Version 7 02-Sep-2003

// Input Error Checking ******************************************************
error(nargchk(1,2,nargin));
if nargin < 2 || isempty(TolInt), TolInt = 0.01*sqrt(eps); end
// End (Input Error Checking) ************************************************

y = abs(x-round(x)) < TolInt;
//--------------------------------------------------------------------------
endfunction

function opt = optstructchk(optdef,optin)
//OPTSTRUCTCHK Validate 'opt' structure.
//    opt = optstructchk(optdef,varargin)
// optdef = default opt structure
//        = MFUN('defaults'), to get default for function MFUN
//  optin = input opt cell from varargin
//    opt = structure, if valid
//        = string, if error
//
// Note: Option values are not vailidated, just its structure
// Used in MINISUMLOC and MCNF

// Copyright (c) 1994-2003 by Michael G. Kay
// Matlog Version 7 02-Sep-2003

if length(optin) == 1, optin = optin{:}; end
opt = optdef;

if iscell(optin)
   if length(optin) == 1 || (length(optin)/2 ~= floor(length(optin)/2))
      opt = 'Cell array ''opt'' must contain field-value pairs.';
   else
      for i = 1:2:length(optin)-1
         if ~isfield(optdef,optin{i})
            opt = 'Field in cell array ''opts'' not valid.';
            return
         elseif ~isempty(optin{i+1})
            opt = setfield(opt,optin{i},optin{i+1});
         end
      end
   end
elseif isstruct(optin)
   fnin = fieldnames(optin); fn = fieldnames(optdef);
   if ~strcmpi(strcat(fnin{:}),strcat(fn{:}))
      opt = '''opt'' not valid options structure.';
   else
      opt = optin;
      for i = 1:length(fnin)
         if isempty(getfield(opt,fnin{i}))
            opt = '''opt'' structure cannot contain an empty field.';
            return
         end
      end
   end
elseif ~isempty(optin)
   opt = '''opt'' must be an options structure or cell array.';
end
//------------------------------------------------------------------------
endfunction
