function [ret_minval,final_xatmin,history] = direct(Problem,bounds,opts,varargin)
// Function   : Direct Version 4.0
// Written by : Dan Finkel (definkel@unity.ncsu.edu)
// Created on : 01/27/2003
// Last Update: 06/21/2004
// Purpose    : Direct optimization algorithm.
//
// This code comes with no guarantee or warranty of any kind
//
//  If unfamiliar with DIRECT, user is recommended to manual
//  that accompanies this code, and can be found at:
//       www4.ncsu.edu/~definkel/research/index.html
//
//
//          [x,fmin,history] = Direct(Problem,bounds,opts,varargin)
//
// Parameters:
//         IN: Problem   - Structure containing problem
//                    Problem.f              = Objective function handle
//
//                    NOTE: If you problem has no constraints (other than
//                          those on the bounds, this is the only field you
//                          need to add to Problem
//
//                    Problem.numconstraints = number of constraints
//                    Problem.constraint(i).func    = i-th constraint handle
//                    Problem.constraint(i).penalty = penalty parameter for
//                                                    i-th constraint
//                    Note: If f returns objective function AND constraints
//                          then set Problem.constraint(1).func = Problem.f
//             bounds    - an n x 2 vector of the lower and upper bounds.
//                         The first column is the lower bounds, and the second
//                         column contains the upper bounds
//             opts      - (optional) MATLAB structure.
//                    opts.ep        = Jones factor                      (default is 1e-4)
//                    opts.maxevals  = max. number of function evals     (default is 2000)
//                    opts.maxits    = max. number of iterations         (default is 1000)
//                    opts.maxdeep   = max. number of rect. divisions    (default is 100)
//                    opts.testflag  = 1 if globalmin known, 0 otherwise (default is 0)
//                    opts.showits   = 1 if disp. stats shown, 0 oth.
//                                      (default is 1)
//                    opts.globalmin = globalmin (if known)
//                                      (default is 0)
//                    opts.tol       = tolerance for term. if tflag=1
//                                      (default is 0.01)
//                    opts.impcons   = turns on implicit constraint capability
//                                      (default is 0)
//                                     If set to one, objective function
//                                     is expected to return a flag which represents
//                                     the feasibility of the point sampled
//             varargin  - (optional) additional arguements to be passed to
//                         objective function
//
//                NOTE: If opts.tflag == 0, maxevals, maxevals and maxdeep are ignored.
//                      DIRECT will stop when the absolute error is less
//                      than tol. Also, preallocation will not occur, and the algorithm
//                      can run slower than if opts.tflag == 1
//                NOTE: opts.maxevals is an approximate stopping condition. DIRECT will
//                      exceed this budget by a slight amount
//
//        OUT: minval    -  minimum value found
//             xatmin    - (optional) location of minimal value
//             history      - (optional) array of iteration historyory, useful for tables and plots
//                The three columns are iteration, fcn evals, and min value found.
//
//                Direct may be called by
//
//                minval = Direct(Problem,bounds);
//                          or
//                with any variation of the optional arguments

//------------------------------------------------------------------//
//
// Implementation taken from:
// D.R. Jones, C.D. Perttunen, and B.E. Stuckman. "Lipschitzian
// Optimization Without the Lipschitz Constant". Journal of
// Optimization Theory and Application, 79(1):157-181, October 1993
//
//------------------------------------------------------------------//
//------------------------------------------------------------------//
Minsize=1e-4;
//------------------------------------------------------------------//
//------------------------------------------------------------------//
//-- Initialize the variables --------------------------------------//
lengths = [];c = [];fc = [];
con = [];szes = [];feas_flags=[];
om_lower     = bounds(:,1);
om_upper     = bounds(:,2);
fcncounter   = 0;
perror       = 0;
itctr        = 1;
done         = 0;
g_nargout    = nargout;
n            = size(bounds,1);

// Determine option values
if nargin<3, opts.tol=0.01; end
if (nargin>=3) && (length(opts)==0), opts.tol=0.01; end

if isfield(opts,'maxits')
  maxits=opts.maxits;
else 
  maxits=1000;          // maximum of iterations
end
if isfield(opts,'maxevals')
  maxevals=opts.maxevals;
else
  maxevals=2000;         // maximum # of function evaluations
end
if isfield(opts,'maxdeep')
  maxdeep=opts.maxdeep;
else
  maxdeep=100;           // maximum number of side divisions
end
if isfield(opts,'testflag')
  testflag=opts.testflag;
else
  testflag=0;          // terminate if within a relative tolerence of f_opt
end
if isfield(opts,'globalmin')
 globalmin=opts.globalmin;
else    
 globalmin=0;          // minimum value of function
end
if isfield(opts,'ep')
 ep=opts.ep;
else    
 ep=1e-4;       // global/local weight parameter.
end
if isfield(opts,'tol')
 tol=opts.tol;
else    
 tol=0.01;       // allowable relative error if f_reach is set
end
if isfield(opts,'showits')
 showits=opts.showits;
else    
 showits=1;          // print iteration stats
end
if isfield(opts,'impcons')
 impcons=opts.impcons;
else    
 impcons=0;          // flag for using implicit constraint handling
end
if isfield(opts,'pert')
 pert=opts.pert;
else    
 pert=1e-6;         // pertubation for implicit constraint handling
end

// 'maxflag',   0, ...         // set to 1 for max problems, 0 for min problems
 //'sizeconst', 0.5,...        // constant on rectangle size function
// 'distance',  1,...          // 1/0 for distance/volume measure of size
// 'minlength', 1e-4,...       // stop if best rectangle has all sides 1ess than this
// 'minevals',  0,...          // but must evaluate at least this many points

theglobalmin = globalmin;
tflag        = testflag;

//-- New 06/08/2004 Pre-allocate memory for storage vectors
if tflag == 0
    lengths    = zeros(n,maxevals + floor(.10*maxevals));
    c          = lengths;
    fc         = zeros(1,maxevals + floor(.10*maxevals));
    szes       = fc;
    con        = fc;
    feas_flags = fc;
end

//-- Call DIRini ---------------------------------------------------//
[thirds , lengths, c , fc, con, feas_flags minval,xatmin,perror,...
        history,szes,fcncounter,calltype] =...
        DIRini(Problem,n,bounds(:,1),bounds(:,2),...
        lengths,c,fc,con, feas_flags, szes,...
        theglobalmin,maxdeep,tflag,g_nargout, impcons, varargin(:));

ret_minval = minval;
ret_xatmin = xatmin;
//-- MAIN LOOP -----------------------------------------------------//
minval = fc(1) + con(1);
while perror > tol
   //-- Create list S of potentially optimal hyper-rectangles
   S = find_po(fc(1:fcncounter)+con(1:fcncounter),...
       lengths(:,1:fcncounter),minval,ep,szes(1:fcncounter));

   //-- Loop through the potentially optimal hrectangles -----------//
   //-- and divide -------------------------------------------------//
   for i = 1:size(S,2)
      [lengths,fc,c,con,feas_flags,szes,fcncounter,success] = ...
          DIRdivide(bounds(:,1),bounds(:,2),Problem,S(1,i),thirds,lengths,...
          fc,c,con,feas_flags,fcncounter,szes,impcons,calltype,varargin(:));
   end

   //-- update minval, xatmin --------------------------------------//
   [minval,fminindex] =  min(fc(1:fcncounter)+con(1:fcncounter));
   penminval = minval + con(fminindex);
   xatmin = (om_upper - om_lower).*c(:,fminindex) + om_lower;
   if (con(fminindex) > 0)|(feas_flags(fminindex) ~= 0)
       //--- new minval is infeasible, don't do anything
   else
       //--- update return values
       ret_minval = minval;
       ret_xatmin = xatmin;
   end

   //--see if we are done ------------------------------------------//
   if tflag == 1
      //-- Calculate error if globalmin known
      if theglobalmin ~= 0
          perror = 100*(minval - theglobalmin)/abs(theglobalmin);
      else
          perror = 100*minval;
      end
   else
      //-- Have we exceeded the maxits?
      if itctr >= maxits
         disp('Exceeded max iterations. Increase maxits')
         done = 1;
      end
      
      // MODIFICAÇÃO : critério de parada : tamanho do quadrado!
      if szes(fminindex)<Minsize
          done = 1;
      end
      //-- Have we exceeded the maxevals?
      if fcncounter > maxevals
         disp('Exceeded max fcn evals. Increase maxevals')
         done = 1;
      end
      if done == 1
         perror = -1;
      end
   end
   if max(max(lengths)) >= maxdeep
      //-- We've exceeded the max depth
      disp('Exceeded Max depth. Increse maxdeep')
      perror = -1;
   end
   if g_nargout == 3
      //-- Store History
      maxhist = size(history,1);
      history(maxhist+1,1) = itctr;
      history(maxhist+1,2) = fcncounter;
      history(maxhist+1,3) = minval;
  end

  //-- New, 06/09/2004
  //-- Call replaceinf if impcons flag is set to 1
  if impcons == 1
      fc = replaceinf(lengths(:,1:fcncounter),c(:,1:fcncounter),...
          fc(1:fcncounter),con(1:fcncounter),...
          feas_flags(1:fcncounter),pert);
  end

  //-- show iteration stats
  if showits == 1
    if  (con(fminindex) > 0) || (feas_flags(fminindex) == 1)
        printf('Iter: %4i   f_min: %15.10f*    fn evals: %8i\n',...
         itctr,minval,fcncounter);
    else
        printf('Iter: %4i   f_min: %15.10f    fn evals: %8i\n',...
         itctr,minval,fcncounter);
    end
  end
  itctr  = itctr + 1;
end

//-- Return values
if g_nargout == 2
    //-- return x*
    final_xatmin = ret_xatmin;
elseif g_nargout == 3
    //-- return x*
    final_xatmin = ret_xatmin;

    //-- chop off 1st row of history
    history(1:size(history,1)-1,:) = history(2:size(history,1),:);
    history = history(1:size(history,1)-1,:);
end
return
endfunction

//------------------------------------------------------------------//
// Function:   DIRini                                               //
// Written by: Dan Finkel                                           //
// Created on: 10/19/2002                                           //
// Purpose   : Initialization of Direct                             //
//             to eliminate storing floating points                 //
//------------------------------------------------------------------//
function [l_thirds,l_lengths,l_c,l_fc,l_con, l_feas_flags, minval,xatmin,perror,...
        history,szes,fcncounter,calltype] = DIRini(Problem,n,a,b,...
        p_lengths,p_c,p_fc,p_con, p_feas_flags, p_szes,theglobalmin,...
        maxdeep,tflag,g_nargout,impcons,varargin)

l_lengths    = p_lengths;
l_c          = p_c;
l_fc         = p_fc;
l_con        = p_con;
l_feas_flags = p_feas_flags;
szes         = p_szes;


//-- start by calculating the thirds array
//-- here we precalculate (1/3)^i which we will use frequently
l_thirds(1) = 1/3;
for i = 2:maxdeep
   l_thirds(i) = (1/3)*l_thirds(i-1);
end

//-- length array will store # of slices in each dimension for
//-- each rectangle. dimension will be rows; each rectangle
//-- will be a column

//-- first rectangle is the whole unit hyperrectangle
l_lengths(:,1) = zeros(n,1);

//01/21/04 HACK
//-- store size of hyperrectangle in vector szes
szes(1,1) = 1;

//-- first element of c is the center of the unit hyperrectangle
l_c(:,1) = ones(n,1)/2;

//-- Determine if there are constraints
calltype = DetermineFcnType(Problem,impcons);

//-- first element of f is going to be the function evaluated
//-- at the center of the unit hyper-rectangle.
//om_point   = abs(b - a).*l_c(:,1)+ a;
//l_fc(1)    = feval(f,om_point,varargin{:});
[l_fc(1),l_con(1), l_feas_flags(1)] = ...
    CallObjFcn(Problem,l_c(:,1),a,b,impcons,calltype,varargin(:));
fcncounter = 1;


//-- initialize minval and xatmin to be center of hyper-rectangle
xatmin = l_c(:,1);
minval   = l_fc(1);
if tflag == 1
    if theglobalmin ~= 0
        perror = 100*(minval - theglobalmin)/abs(theglobalmin);
    else
        perror = 100*minval;
    end
else
   perror = 2;
end

//-- initialize history
//if g_nargout == 3
    history(1,1) = 0;
    history(1,2) = 0;
    history(1,3) = 0;
//end
endfunction

//------------------------------------------------------------------//
// Function   :  find_po                                            //
// Written by :  Dan Finkel                                         //
// Created on :  10/19/2002                                         //
// Purpose    :  Return list of PO hyperrectangles                  //
//------------------------------------------------------------------//
function rects = find_po(fc,lengths,minval,ep,szes)

//-- 1. Find all rects on hub
diff_szes = sum(lengths,1);
tmp_max = max(diff_szes);
j=1;
sum_lengths = sum(lengths,1);
for i =1:tmp_max+1
    tmp_idx = find(sum_lengths==i-1);
    [tmp_n, hullidx] = min(fc(tmp_idx));
    if length(hullidx) > 0
        hull(1,j) = tmp_idx(hullidx);
        j=j+1;
        //-- 1.5 Check for ties
        ties = find(abs(fc(tmp_idx)-tmp_n) <= 1e-13);
        if length(ties) > 1
            mod_ties = find(tmp_idx(ties) ~= hull(j-1));
            hull = [hull tmp_idx(ties(mod_ties))];
            j = length(hull)+1;
        end
    end
end
//-- 2. Compute lb and ub for rects on hub
lbound = calc_lbound(lengths,fc,hull,szes)';
ubound = calc_ubound(lengths,fc,hull,szes)';

//-- 3. Find indeces of hull who satisfy
//--    1st condition
maybe_po = find(lbound-ubound <= 0);

//-- 4. Find indeces of hull who satisfy
//--    2nd condition
t_len  = length(hull(maybe_po));
if minval ~= 0
    po = find((minval-fc(hull(maybe_po)))./abs(minval) +...
        szes(hull(maybe_po)).*ubound(maybe_po)./abs(minval) >= ep);
else
    po = find(fc(hull(maybe_po)) -...
        szes(hull(maybe_po)).*ubound(maybe_po) <= 0);
end
final_pos      = hull(maybe_po(po));

rects = [final_pos; szes(final_pos)];
return
endfunction

//------------------------------------------------------------------//
// Function   :  calc_ubound                                        //
// Written by :  Dan Finkel                                         //
// Created on :  10/19/2002                                         //
// Purpose    :  calculate the ubound used in determing potentially //
//               optimal hrectangles                                //
//------------------------------------------------------------------//
function ub = calc_ubound(lengths,fc,hull,szes)

hull_length  = length(hull);
hull_lengths = lengths(:,hull);
for i =1:hull_length
    tmp_rects = find(sum(hull_lengths,1)<sum(lengths(:,hull(i))));
    if length(tmp_rects) > 0
        tmp_f     = fc(hull(tmp_rects));
        tmp_szes  = szes(hull(tmp_rects));
        tmp_ubs   = (tmp_f-fc(hull(i)))./(tmp_szes-szes(hull(i)));
        ub(i)        = min(tmp_ubs);
    else
        ub(i)=1.976e14;
    end
end
return
endfunction

//------------------------------------------------------------------//
// Function   :  calc_lbound                                        //
// Written by :  Dan Finkel                                         //
// Created on :  10/19/2002                                         //
// Purpose    :  calculate the lbound used in determing potentially //
//               optimal hrectangles                                //
//------------------------------------------------------------------//
function lb = calc_lbound(lengths,fc,hull,szes)

hull_length  = length(hull);
hull_lengths = lengths(:,hull);
for i = 1:hull_length
    tmp_rects = find(sum(hull_lengths,1)>sum(lengths(:,hull(i))));
    if length(tmp_rects) > 0
        tmp_f     = fc(hull(tmp_rects));
        tmp_szes  = szes(hull(tmp_rects));
        tmp_lbs   = (fc(hull(i))-tmp_f)./(szes(hull(i))-tmp_szes);
        lb(i)     = max(tmp_lbs);
    else
        lb(i)     = -1.976e14;



    end
end
return
endfunction

//------------------------------------------------------------------//
// Function   :  DIRdivide                                          //
// Written by :  Dan Finkel                                         //
// Created on :  10/19/2002                                         //
// Purpose    :  Divides rectangle i that is passed in              //
//------------------------------------------------------------------//
function [lengths,fc,c,con,feas_flags,szes,fcncounter,pass] = ...
    DIRdivide(a,b,Problem,index,thirds,p_lengths,p_fc,p_c,p_con,...
    p_feas_flags,p_fcncounter,p_szes,impcons,calltype,varargin)

lengths    = p_lengths;
fc         = p_fc;
c          = p_c;
szes       = p_szes;
fcncounter = p_fcncounter;
con        = p_con;
feas_flags = p_feas_flags;

//-- 1. Determine which sides are the largest
li     = lengths(:,index);
biggy  = min(li);
ls     = find(li==biggy)';
lssize = length(ls);
j = 0;

//-- 2. Evaluate function in directions of biggest size
//--    to determine which direction to make divisions
oldc       = c(:,index);
delta      = thirds(biggy+1);
newc_left  = oldc(:,ones(1,lssize));
newc_right = oldc(:,ones(1,lssize));
f_left     = zeros(1,lssize);
f_right    = zeros(1,lssize);
for i = 1:lssize
    lsi               = ls(i);
    newc_left(lsi,i)  = newc_left(lsi,i) - delta;
    newc_right(lsi,i) = newc_right(lsi,i) + delta;
    [f_left(i), con_left(i), fflag_left(i)]    = CallObjFcn(Problem,newc_left(:,i),a,b,impcons,calltype,varargin(:));
    [f_right(i), con_right(i), fflag_right(i)] = CallObjFcn(Problem,newc_right(:,i),a,b,impcons,calltype,varargin(:));
    fcncounter = fcncounter + 2;
end
w = [min(f_left, f_right)' ls];

//-- 3. Sort w for division order
[V,order] = gsort(w,"r","i");

//-- 4. Make divisions in order specified by order
for i = 1:size(order,1)

   newleftindex  = p_fcncounter+2*(i-1)+1;
   newrightindex = p_fcncounter+2*(i-1)+2;
   //-- 4.1 create new rectangles identical to the old one
   oldrect = lengths(:,index);
   lengths(:,newleftindex)   = oldrect;
   lengths(:,newrightindex)  = oldrect;

   //-- old, and new rectangles have been sliced in order(i) direction
   lengths(ls(order(i,1)),newleftindex)  = lengths(ls(order(i,1)),index) + 1;
   lengths(ls(order(i,1)),newrightindex) = lengths(ls(order(i,1)),index) + 1;
   lengths(ls(order(i,1)),index)         = lengths(ls(order(i,1)),index) + 1;

   //-- add new columns to c
   c(:,newleftindex)  = newc_left(:,order(i));
   c(:,newrightindex) = newc_right(:,order(i));

   //-- add new values to fc
   fc(newleftindex)  = f_left(order(i));
   fc(newrightindex) = f_right(order(i));

   //-- add new values to con
   con(newleftindex)  = con_left(order(i));
   con(newrightindex) = con_right(order(i));

   //-- add new flag values to feas_flags
   feas_flags(newleftindex)  = fflag_left(order(i));
   feas_flags(newrightindex) = fflag_right(order(i));

   //-- 01/21/04 Dan Hack
   //-- store sizes of each rectangle
   szes(1,newleftindex)  = 1/2*norm((1/3*ones(size(lengths,1),1)).^(lengths(:,newleftindex)));
   szes(1,newrightindex) = 1/2*norm((1/3*ones(size(lengths,1),1)).^(lengths(:,newrightindex)));
end
szes(index) = 1/2*norm((1/3*ones(size(lengths,1),1)).^(lengths(:,index)));
pass = 1;

return
endfunction

//------------------------------------------------------------------//
// Function   :  CallConstraints                                    //
// Written by :  Dan Finkel                                         //
// Created on :  06/07/2004                                         //
// Purpose    :  Evaluate Constraints at pointed specified          //
//------------------------------------------------------------------//
function ret_value = CallConstraints(Problem,x,a,b,varargin)

//-- Scale variable back to original space
point = abs(b - a).*x+ a;

ret_value = 0;
if isfield(Problem,'constraint')
    if ~isempty(Problem.constraint)
        for i = 1:Problem.numconstraints
            if length(Problem.constraint(i).func) == length(Problem.f)
                if double(Problem.constraint(i).func) == double(Problem.f)                    
                    //-- Dont call constraint; value was returned in obj fcn
                    con_value = 0;
                else
                    con_value = feval(Problem.constraint(i).func,point,varargin(:));
                end
            else
                con_value = feval(Problem.constraint(i).func,point,varargin(:));
            end
            if con_value > 0
                //-- Infeasible, punish with associated pen. param
                ret_value = ret_value + con_value*Problem.constraint(i).penalty;
            end
        end
    end
end
return
endfunction

//------------------------------------------------------------------//
// Function   :  CallObjFcn                                         //
// Written by :  Dan Finkel                                         //
// Created on :  06/07/2004                                         //
// Purpose    :  Evaluate ObjFcn at pointed specified               //
//------------------------------------------------------------------//
function [fcn_value, con_value, feas_flag] = ...
    CallObjFcn(Problem,x,a,b,impcon,calltype,varargin)

con_value = 0;
feas_flag = 0;

//-- Scale variable back to original space
point = abs(b - a).*x+ a;

if calltype == 1
    //-- No constraints at all
    fcn_value = evstr(Problem.f+"(point)",varargin(:));
elseif calltype == 2
    //-- f returns all constraints
    [fcn_value, cons] = evstr(Problem.f+"(point)",varargin(:));
    for i = 1:length(cons)
        if cons > 0
            con_value = con_value + Problem.constraint(i).penalty*cons(i);
        end
    end
elseif calltype == 3
    //-- f returns no constraint values
    fcn_value = evstr(Problem.f+"(point)",varargin(:));
    con_value = CallConstraints(Problem,x,a,b,varargin(:));
elseif calltype == 4
    //-- f returns feas flag
    [fcn_value,feas_flag] = evstr(Problem.f+"(point)",varargin(:));
elseif calltype == 5
    //-- f returns feas flags, and there are constraints
    [fcn_value,feas_flag] = evstr(Problem.f+"(point)",varargin(:));
    con_value = CallConstraints(Problem,x,a,b,varargin(:));
end
if feas_flag == 1
    fcn_value = 10^9;
    con_value = 0;
end
return
endfunction

//------------------------------------------------------------------//
// Function   :  replaceinf                                         //
// Written by :  Dan Finkel                                         //
// Created on :  06/09/2004                                         //
// Purpose    :  Assign R. Carter value to given point              //
//------------------------------------------------------------------//
function fcn_values = replaceinf(lengths,c,fc,con,flags,pert)

//-- Initialize fcn_values to original values
fcn_values = fc;

//-- Find the infeasible points
infeas_points = find(flags == 1);

//-- Find the feasible points
feas_points   = find(flags == 0);

//-- Calculate the max. value found so far
if ~isempty(feas_points)
    maxfc = max(fc(feas_points) + con(feas_points));
else
    maxfc = max(fc + con);
end

for i = 1:length(infeas_points)
    if isempty(feas_points)
        //-- no feasible points found yet
        found_points = [];found_pointsf = [];
        index = infeas_points(i);
    else
        index = infeas_points(i);

        //-- Initialize found points to be entire set
        found_points  = c(:,feas_points);
        found_pointsf = fc(feas_points) + con(feas_points);

        //-- Loop through each dimension, and find points who are close enough
        for j = 1:size(lengths,1)
            neighbors = find(abs(found_points(j,:) - c(j,index)) <= ...
                3^(-lengths(j,index)));
            if ~isempty(neighbors)
                found_points  = found_points(:,neighbors);
                found_pointsf = found_pointsf(neighbors);
            else
                found_points = [];found_pointsf = [];
                break;
            end
        end
    end

    //-- Assign Carter value to the point
    if ~isempty(found_pointsf)
        //-- assign to index the min. value found + a little bit more
        fstar = min(found_pointsf);
        if fstar ~= 0
            fcn_values(index) = fstar + pert*abs(fstar);
        else
            fcn_values(index) = fstar + pert*1;
        end
    else
        fcn_values(index) = maxfc+1;
        maxfc             = maxfc+1;
    end
end
return
endfunction

//------------------------------------------------------------------//
// Function   :  DetermineFcnType                                   //
// Written by :  Dan Finkel                                         //
// Created on :  06/25/2004                                         //
// Purpose    :  Determine how constraints are handled              //
//------------------------------------------------------------------//
function retval = DetermineFcnType(Problem,impcons)

retval = 0;
if (~isfield(Problem,'constraint'))&(~impcons)
    //-- No constraints at all
    retval = 1;
end
if isfield(Problem,'constraint')
    //-- There are explicit constraints. Next determine where
    //-- they are called
    if ~isempty(Problem.constraint)
        if length(Problem.constraint(1).func) == length(Problem.f)
            //-- Constraint values may be returned from objective
            //-- function. Investigate further
            if double(Problem.constraint(1).func) == double(Problem.f)
                //-- f returns constraint values
                retval = 2;
            else
                //-- f does not return constraint values
                retval = 3;
            end
        else
            //-- f does not return constraint values
            retval = 3;
        end
    else
        if impcons
            retval = 0;
        else
            retval = 1;
        end
    end
end

if (impcons)
    if ~retval
        //-- only implicit constraints
        retval = 4;
    else
        //-- both types of constraints
        retval = 5;
    end
end
endfunction

//------------------------------------------------------------------//
// GETOPTS Returns options values in an options structure
// USAGE
//   [value1,value2,...]=getopts(options,field1,default1,field2,default2,...)
// INPUTS
//   options  : a structure variable
//   field    : a field name
//   default  : a default value
// OUTPUTS
//   value    : value in the options field (if it exists) or the default value
//
// Variables with the field names will be created in the caller's workspace
// and set to the value in the option variables field (if it exists) or to the
// default value.
//
// Example called from a function:
//   getopts(options,'tol',1e-8,'maxits',100);
// where options contains the single field 'tol' with value equal to 1
// The function have two variable defined in the local workspace, tol with a
// value of 1 and maxits with a value of 100.
//
// If options contains a field name not in the list passed to getopts, a
// warning is issued.
//
//
// Many thanks to the author of this function,
// Paul Fackler (pfackler@ncsu.edu)
//
//
//------------------------------------------------------------------//
function varargout=getopts(options,varargin)
K=fix(nargin/2);
if nargin/2==K
  error('fields and default values must come in pairs')
end
if typeof(options) == "st", optstruct=1; else, optstruct=0; end
varargout=cell(K,1);
k=0;
ii=1;
for i=1:K
  if optstruct && isfield(options,varargin(ii))
    assignin('caller',varargin{ii},getfield(options,varargin(ii)));
    k=k+1;
  else
    assignin('caller',varargin(ii),varargin(ii+1));
  end
  ii=ii+2;
end

if optstruct && k~=size(fieldnames(options),1)
  warning('options variable contains improper fields')
end

return
//------------------------------------------------------------------//
// Versions  : 1.0 - 1st successful implemenation of DIRect
//           : 2.0 - Removed floating point arithmetic
//                   duplicated Table 5 of Jones et al.
//           : 2.1 - increased speed by storing size calcs.
//           : 2.2 - utitilized linked lists to increase speed
//           : 2.3 - rewrote ubound to increase speed
//           : 2.4 - rewrote lbound to increase speed
//           : 2.5 - removed call to calcsize
//           : 2.6 - added check_for_ties
//           : 2.7 - rewrote check_for_ties to compare fp correctly
//           : 2.8 - changed output arguments, rewrote help
//           : 3.0 - simplified input/output. Put on web.
//           : 3.1 - Performanced Tuned! Tremendous speed increase
//           : 3.2 - Removed llists; performance tuned
//                   Many thanks to Ray Muzic and Paul Fackler
//                   for their suggestions to improve this code
//           : 4.0 - Sped up code, and added 2 constraint handling
//                   mechanisms.
//------------------------------------------------------------------//
endfunction
