function [wo, nNLPs] = minlp_nm(NLP,n,m,nh,fig,x0,y0,xlb,xub,M,epsR,epsI,epsC,epsZ)
%     * * * Solution of Mixed Integer Non Linear Programming * * *
%    Z = min(x,y) F(x,y)
%    s.t.         G(x,y) <= 0
%                 x in X in R^n
%                 y in Y = {0,1}^m
%
%  usage: [wo, nIter] = minlp_nm(NLP,n,m,nh,fig,x0,y0,xlb,xub,M,epsR,epsI,epsC,epsZ)
%
%  wo      : best solution found [x,y,u,Z]
%            where u is the maximum constraint violation
%  St      : matrix = (number of constraints, number of LPs) for each NLP
%  NLP     : M file describing the NLP: [F,G] = nlp(x,y)
%  n       : number of continuous variables
%  m       : number of binary variables
%  nh      : number of equality constraints
%  fig     : = 1 : plot constraints surfaces. = 0 : no plots (default).
%  x0      : initial guess for x (default = zeros(1,n)). Returns first solution.
%  y0      : initial guess for y (default = round(rand(1,m)))
%  xlb     : lower bound for x (default = zeros(1,n))
%  xub     : upper bound for x (default = inf * ones(1,n))
%  M > 0   : constraint violation penalization (default = 1e3)
%  epsR    : real variable tolerance (default = 1e-4)
%  epsI    : integer variable tolerance (default = 1e-4)
%  epsC    : constraint violation tolerance (default = 1e-6)
%  epsZ    : objective function tolerance (default = 1e-4)

% adaptado de minlpn.m

if nargin < 3
  disp('minimal usage: minlp(NLP,n,m)');
  return;
end
if nargin < 4
   nh = 0;
end
if (nh < 0)
   disp('MINLPN Error: nh cannot be lower than zero');
   return;
end
if nargin < 5
  fig = 0;
end
if nargin < 6
  x0 = zeros(1,n);
end
if nargin < 7
  y0 = round(rand(1,m));
end
if nargin < 8
  xlb = zeros(1,n);
end
if nargin < 9
  xub = inf * ones(1,n);
end
if nargin < 10
  M = 1e3;
end
if nargin < 11
  epsR = 1e-4;
end
if nargin < 12
  epsI = 1e-4;
end
if nargin < 13
  epsC = 1e-6;
end
if nargin < 14
  epsZ = 1e-4;
end

if isempty(nh)
   nh = 0;
end
if isempty(fig)
  fig = 0;
end
if isempty(x0)
  x0 = zeros(1,n);
end
if isempty(y0)
  y0 = round(rand(1,m));
end
if isempty(xlb)
  xlb = zeros(1,n);
end
if isempty(xub)
  xub = inf * ones(1,n);
end
if isempty(M) | M <= 0
  M = 1e3;
end
if isempty(epsR)
  epsR = 1e-4;
end
if isempty(epsI)
  epsI = 1e-4;
end
if isempty(epsC)
  epsC = 1e-6;
end
if isempty(epsZ)
  epsZ = 1e-4;
end

ylb = zeros (1,m);
yub = ones(1,m);

op = foptions;
op(1) = 0;        % 1 - display constr intermediary results; 0 - nothing
op(2) = epsR;
op(3) = epsZ;
op(4) = epsC;
op(13) = nh;
op(14) = n*100;  % max function evaluation for complex

unFeasible = 1; % flag para marcar se o problema inicial é viável ou não
nNLPs = 0; %número de NLPs resolvidos
maxInit = 10; %número máximo de tentativas de inicializar
flag = 0;
zup = inf;
no = 1;
k = no;

while unFeasible
   
   unFeasible = 0;
   
   if(nNLPs >0) % trocando o chute inicial
      x0 = rand(1,n);
      y0 = round(rand(1,m));
   end
   
   % NLP with relaxed integrality
   wo=[x0 y0 -inf inf];
   wlb=[xlb ylb];
   wub=[xub yub];
   
   W = []; Tr = [];
   [flag,wo,W,Tr,no,unFeasible] = nlp_nm_internal(NLP,n,m,nh,fig,wo,M,wlb,wub,W,Tr,0,no,epsR,epsI,epsC,epsZ);

   y0 = round(wo(n+1:n+m));
   [dym,j] = max(abs(wo(n+1:n+m)-y0)); %verificando se a solução é inteira em y
   x0 = wo(1:n);
   
   if dym > epsI %se a solução do prob. relaxado é não inteira flag == 1
      flag = 1;
      wo(n+1:n+m) = y0;
      [x0,ot,nS] = complexb('fun_nm',x0,0,xlb,xub,-1,op(3),op(14),[],[],NLP,y0,n,nh,M);
      op(10)=nS;      
      [F, G] = eval([NLP '(x0,y0)']);    
      wo(n+m+1) = max([abs(G(1:nh)) G(nh+1:end) 0]);
      if wo(n+m+1)>0
         wo(n+m+2) = inf;
      else
         wo(n+m+2) = F;
      end
   else
      no = 0; %chegou na solução do problema
   end
   
      
   nNLPs = nNLPs+ 1;
   if(nNLPs > maxInit)
      disp('*******************');
      disp('****MINLPN error****: no feasible solution found for the totally relaxed problem, try increase maxInit.');
      disp('*******************');
      no =0;
      break;
   end
   
end % while (unfeasible)

%r = length(g); % number of constraints

%wo = [x0,y0,u,zu];  % best solution at each iteration

Tr = [0, 1, 0, 0, 1:m];  % binary tree [parent, child, mj, nfix, fix, ~fix]
W = wo;                  % solution of the LPs
Tr(1,3) = j;    % less integer yj

while no
   if flag             % non integer solution
      w = W(no,:);      % create two child problems deleting the parent
      j = Tr(no,3);
      pai = Tr(no,2);
      nfx = Tr(no,4);
      fx = [Tr(no,5:4+nfx),j];
      vr = Tr(no,5+nfx:4+m);
      vr = vr(find(vr-j));
      Tr(no,:) = [pai,pai+1,0,nfx+1,fx,vr];
      no = no + 1;
      W(no,:) = w;
      Tr(no,:) = Tr(no-1,:);
      Tr(no,2) = pai+2;
      j = j + n;
      yj = round(w(j));
      W(no-1,j) = yj;
      W(no,j) = 1 - yj;
   end
   
   % round trip for version 4.2
   % version 5.2:   if no > 1 & Tr(no,1) == Tr(no-1,1)
   
   if no > 1
      xx = Tr(no-1,1);
   else
      xx = -1;
   end
   
   if Tr(no,1) == xx     % last two child problems with same parent
      old_no = no;
      [flag,wo,W,Tr,no,unFeasible] = nlp_nm_internal(NLP,n,m,nh,fig,wo,M,wlb,wub,W,Tr,no,no,epsR,epsI,epsC,epsZ);
      nNLPs = nNLPs+1;
      if flag
         k = old_no - 1;
         [flag,wo,W,Tr,no,unFeasible] = nlp_nm_internal(NLP,n,m,nh,fig,wo,M,wlb,wub,W,Tr,k,no,epsR,epsI,epsC,epsZ);
         nNLPs = nNLPs+1;
         if flag & old_no == no
            if W(old_no,n+m+2) > W(k,n+m+2)
               W([k old_no],:) = W([old_no k],:); % switch last two nodes
               Tr([k old_no],:) = Tr([old_no k],:);
            end
         end
      end
   else
      [flag,wo,W,Tr,no,unFeasible] = nlp_nm_internal(NLP,n,m,nh,fig,wo,M,wlb,wub,W,Tr,no,no,epsR,epsI,epsC,epsZ);
      nNLPs = nNLPs+1;
   end
end


%% testando se foi encontrada solução
if (unFeasible)
   wo(n+m+2) = inf;
   disp('*******************');
   disp('****MINLPN error****: NO feasible solution found.');
   disp('*******************');
else
   disp('****MINLPN****: solution found.');
end
