function [wo,St,x0,y0] = minlp(NLP,GRAD,n,m,rx,fig,x0,y0,xlb,xub,M,epsR,epsI,epsC,epsZ)
%     * * * Solution of MINLP * * *
%   Z = min(y,x) F(x,y) = C'y + f(x)
%    s.t. G(x,y) = By + g(x) <= 0
%         x in X in R^n
%         y in Y = {0,1}^m
%
%  usage: [wo,St,x0,y0] = minlp(NLP,GRAD,n,m,rx,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)
%            with the first line: if size(x,2) > n,  y = x(n+1:n+m); end
%  GRAD    : M file for gradients: [Fx,Gx,C,B] = grad(x,y)
%            with the first line: if size(x,2) > n,  y = x(n+1:n+m); end
%            and last lines: if size(x,2) > n, Fx = [Fx C]; Gx = [Gx B]; end
%                            Gx = Gx';
%  rx      : = 1 : relaxed initial NLP. = 0 : unrelaxed initial NLP (default)
%  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)
%

% authors: Argimiro R. Secchi and Elaine C. Pereira
%          Chemical Engineering Department
%          Applied Math Department
%          Universidade Federal do Rio Grande do Sul, Brasil
%          arge@enq.ufrgs.br
%
% version 1.0: 15/mar/1999

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

if isempty(rx)
  rx = 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);

%   * * * Solution of initial NLP * * *
%        Z = min(x) C'y0 + f(x) + Mu
%        s.t. By0 + g(x) <= u
%             x in X
%             y0 fixed
%             u in R+

if rx           % NLP with relaxed integrality
  x0=[x0 y0];
  wlb=[xlb ylb];
  wub=[xub yub];
else
  wlb=xlb;
  wub=xub;
end

op = foptions;
op(1) = 1;        % 1 - display some results; 0 - nothing
op(2) = epsR;
op(3) = epsZ;
op(4) = epsC;
op(13) = 0;
op(14) = 150;  % max function evaluation

% constr returns the best solution, even violating the constraints
% then u does not need be computed explicitly.

[x0,op,lambda] = constr(NLP,x0,op,wlb,wub,GRAD,y0);

if rx
  y0 = round(x0(n+1:n+m));
  dym = max(abs(x0(n+1:n+m)-y0));
  x0 = x0(1:n);
  if dym > epsI
    [x0,op,lambda] = constr(NLP,x0,op,xlb,xub,GRAD,y0);
  end   
end

% getting the values of f, g, and their gradients
[f_nlp,g_nlp] = eval([NLP '(x0,y0)']);
u = max(g_nlp);             % u > 0 means constraint violation
[gf_nlp,gg_nlp,C,B] = eval([GRAD '(x0,y0)']);
f = f_nlp - C * y0';
g = g_nlp' - B * y0';
gf = gf_nlp;
gg = gg_nlp';

if u > epsC                 % Z upper bound
  zu = inf;
else
  zu = f_nlp;
end

r = size(B,1);             % number of constraints
St = [r+1,0];              % matrix of NLPs: number of constraints x LPs

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

A = [gg,B,-ones(r,1),zeros(r,1); gf,C,M,-1];
b = [gg*x0'- g; gf*x0'-f]; 

wlb = [xlb,ylb,0,-inf];
wub = [xub,yub,inf,inf];

if fig, planos(n,A,b,y0,u,xlb,xub,1); end

Tr = [0, 1, 0, 0, 1:m];  % binary tree [parent, child, mj, nfix, fix, ~fix]
W = wo;                  % solution of the LPs
no = 1;                  % number of nodes in the tree

[flag,wo,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,1,no,epsR,epsI,epsC,epsZ,St);

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,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,no,no,epsR,epsI,epsC,epsZ,St);
    if flag
      k = old_no - 1;
      [flag,wo,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,k,no,epsR,epsI,epsC,epsZ,St);
       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,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,no,no,epsR,epsI,epsC,epsZ,St);
  end
end

