function [x,fval,exitflag,output,lambda,nlp,pos]=milprog(f,n,m,A,B,Aeq,Beq,lb,ub,x0,options)
% 
% %MILPROG     Linear programming by exaustive search.                 
% %   X=MILPROG(f,n,m,A,b) solves the mixed integer linear programming problem:
% %        
% %            min f'*X    subject to:   A*X <= b 
% %            x,y
% %            x real
% %            y integer
% %            X=[x;y]
% %            n=number of real variables
% %            m=number of integer variables
% %   X=MILPPROG(f,n,m,A,b,Aeq,beq) solves the problem above while additionally
% %   satisfying the equality constraints Aeq*[x;y] = beq.
% %   
% %   X=MILPROG(f,n,m,A,b,Aeq,beq,LB,UB) defines a set of lower and upper
% %   bounds on the design variables, X, so that the solution is in
% %   the range LB <= X <= UB.  Use empty matrices for LB and UB
% %   if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; 
% %   set UB(i) = Inf if X(i) is unbounded above.
% %
% %   X=MILPROG(f,n,m,A,b,Aeq,beq,LB,UB,X0) sets the starting point to X0.  This
% %   option is only available with the active-set algorithm.  The default
% %   interior point algorithm will ignore any non-empty starting point.
% %
% %   X=MILPROG(f,n,m,A,b,Aeq,Beq,LB,UB,X0,OPTIONS) minimizes with the default 
% %   optimization parameters replaced by values in the structure OPTIONS, an 
% %   argument created with the OPTIMSET function.  See OPTIMSET for details.  
% %   Use options are Display, Diagnostics, TolFun, LargeScale, MaxIter. 
% %   Currently, only 'final' and 'off' are valid values for the parameter 
% %   Display when LargeScale is 'off' ('iter' is valid when LargeScale is 'on').
% %
% %   [X,FVAL]=MILPROG(f,n,m,A,b) returns the value of the objective function at X:
% %   FVAL = f'*X.
% %
% %   [X,FVAL,EXITFLAG] = MILPROG(f,A,b) returns EXITFLAG that 
% %   describes the exit condition of LINPROG.
% %   If EXITFLAG is:
% %      > 0 then LINPROG converged with a solution X.
% %      0   then LINPROG reached the maximum number of iterations without converging.
% %      < 0 then the problem was infeasible or LINPROG failed.
% %
% %   [X,FVAL,EXITFLAG,OUTPUT] = MILPROG(f,n,m,A,b) returns a structure
% %   OUTPUT with the number of iterations taken in OUTPUT.iterations, the type
% %   of algorithm used in OUTPUT.algorithm, the number of conjugate gradient
% %   iterations (if used) in OUTPUT.cgiterations.
% %
% %   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,NLP,POS]=MILPROG(f,n,m,A,b) returns the set of 
% %   Lagrangian multipliers LAMBDA, at the solution: LAMBDA.ineqlin for the 
% %   linear inequalities A, LAMBDA.eqlin for the linear equalities Aeq, 
% %   LAMBDA.lower for LB, and LAMBDA.upper for UB.NLP returns the number of
% %   LINPTOG evaluations.POs returns the position of the the solution, use to
%  %  analize the results.
% %   
% %   NOTE: the LargeScale (the default) version of MILPROG uses LINPROG that use a primal-dual
% %         method. Both the primal problem and the dual problem must be feasible 
% %         for convergence. Infeasibility messages of either the primal or dual, 
% %         or both, are given as appropriate.  The primal problem in standard 
% %         form is 
% %              min f'*x such that A*x = b, x >= 0.
% %         The dual problem is
% %              max b'*y such that A'*y + s = f, s >= 0.
% % Authors: Argimiro R. Secchi &
% %          Marcelo Escobar
% %          Chemical Engineering Department
% %          Applied Math Department
% %          Universidade Federal do Rio Grande do Sul, Brasil
% %          arge@enq.ufrgs.br | escobar@enq.ufrgs.br
% %
% % version 1.0: 15/agosto/2004
% 
% If just 'defaults' passed in, return the default options in X

% If just 'defaults' passed in, return the default options in X

defaultopt = struct('Display','final',...
   'TolFun',1e-8,'Diagnostics','off',...
   'LargeScale','on','MaxIter',85);

if nargin==1 & nargout <= 1 & isequal(f,'defaults')
   x = defaultopt;
   return
end

% Handle missing arguments

 if nargin < 11, options = [];
   if nargin < 10, x0 = []; 
      if nargin < 9, ub = []; 
         if nargin < 8, lb = []; 
            if nargin < 7, Beq = [];
               if nargin < 6, Aeq = [];
               end, end, end, end, end, end
if nargout > 4
   computeLambda = 1;
else 
   computeLambda = 0;
end

%M-Matriz de Possibilidades:
nv=length(f);
ne=size(Aeq);
ni=size(A);
ny=m;
N=[];
M=zeros(2^ny,ny);
for i=1:ny
   eq=[zeros(2^(ny-i),1);ones(2^(ny-i),1)];
   cp=length(eq);
      if cp==2^ny
          N=eq;
      else N=[];
          while cp~= 2^ny
             N=[N;eq];
             cp=length(N);
             end
      end
 M(:,i)=N;
 end
 tol=1e-6;
 if nv==m % Apenas variaveis Binarias
     fv=[];
     pos=[];
    if ne(1)==0 %Apenas Equacoes de Desigualdade
         for i=1:2^ny
         flagv(i)=  sum((A*M(i,:)'-B)>=0);   %flagv=1 implica em violacao de restricao
         if flagv(i)==0
             fv=[fv;f*M(i,:)'];
             pos=[pos;i];
         end
         end
     fval=min(fv);
     x= M(pos(find(fv==min(fv))),:)';
     pos=pos(find(fv==min(fv)));
    else %Equacoes de Igauladdae e Desigualdade
         for i=1:2^ny
         flagv(i)=sum(abs(Aeq*M(i,:)'-Beq)>=tol)+sum((A*M(i,:)'-B)>=0);  %flagv=1 implica em violacao de restricao
         if flagv(i)==0
             fv=[fv;f*M(i,:)'];
             pos=[pos;i];
         end
         end
     fval=min(fv);
     x= M(pos(find(fv==min(fv))),:)';
     pos=pos(find(fv==min(fv)));
    end
     exitflag=1;
     output=[];
     lambda=[];
     nlp=0;
else %Variaveis Binarias e Reais
    if ne(1)==0 %Apenas Equacoes de Desigualdade
        c=f(1:nv-m);
        Am=A(:,1:nv-m);
        lb=lb(1:nv-m);
        ub=ub(1:nv-m);
        op=options;
        fv=[];pos=[];Mxy=[];nlp=0;
        for i=1:2^ny
            Bm=B-f(nv-m+1:end)*M(i,:)';
            [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(c,Am,Bm,Aeq,Beq,lb,ub,x0,op);
            exitflag{i}=EXITFLAG;
            output{i}=OUTPUT;
            lambda{i}=LAMBDA;
            nlp=nlp+(EXITFLAG>0);
            xmy=[X;M(i,:)'];
            Mxy=[Mxy;xmy'];
            FVALy(i)=FVAL+f(nv-m+1:end)*M(i,:)';
            flagv(i)=  sum((A*xmy-B)>=0);
            if flagv(i)==0
               fv=[fv;f*xmy];
               pos=[pos;i];
            end
        end
        fval=min(fv);
        x= Mxy(pos(find(fv==min(fv))),:)';
        pos=pos(find(fv==min(fv)));
        
    else %Equacoes de Igauladdae e Desigualdade 
     Aeqm=[Aeq(:,1:nv-m)  eye(ne(1)) -eye(ne(1))];
     Am=[ A(:,1:nv-m) zeros(ni(1),2*ne(1))];
     lb=[lb(1:nv-m) zeros(1,2*ne(1))];
     ub=[ub(1:nv-m) Inf*ones(1,2*ne(1))];
     c=[f(1:nv-m) ones(1,2*ne(1))];
     op=options;
     fv=[];pos=[];Mxy=[];nlp=0;
     for i=1:2^ny
            Bm=B-f(nv-m+1:end)*M(i,:)';
            Beqm=Beq-f(nv-m+1:end)*M(i,:)';
            [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(c,Am,Bm,Aeqm,Beqm,lb,ub,x0,op);
            exitflag{i}=EXITFLAG;
            output{i}=OUTPUT;
            lambda{i}=LAMBDA;
            nlp=nlp+(EXITFLAG>0);
            xmy=[X;M(i,:)'];
            Mxy=[Mxy;xmy'];
            FV(i)=FVAL;
            FVALy(i)=c(1:nv-m)*X(1:nv-m)+f(nv-m+1:end)*M(i,:)';
            flagv(i)= sum((A*[xmy(1:nv-m);M(i,:)'] -B)>=0) + sum((Aeq*[xmy(1:nv-m);M(i,:)'] -Beq)>=tol);
            if EXITFLAG==1 
                if flagv(i)==0
                   fv=[fv;f*[xmy(1:nv-m);M(i,:)']];
                   pos=[pos;i];
                end
            end
        end
        fval=min(fv);
        MXY=[Mxy(:,1:nv-m) M];
        x= MXY(pos(find(fv==min(fv))),:)';
        pos=pos(find(fv==min(fv)));
    end
end
disp('====================================');
disp('====================================');
disp('RESULTADOS OBTIDOS:')