function [x_opt,f_opt,g_opt,cvg_hst]=ORSopt(func,x_init,x_lb,x_ub,options,consts)
% [x_opt,f_opt,g_opt,cvg_hst]=ORSopt(func,x_init,x_lb,x_ub,options,consts)
%
% ORSopt: Optimized Step Size Randomized Search
% Nonlienar optimization with inequality constraints using Random Search
%
%       minimizes f(x) such that g(x)<0 and x_lb <= x_opt <= x_ub.
%       f is a scalar objective function, x is a vector of design variables, and
%       g is a vector of constraints.
%
% INPUT
% ======
%  func    :   the name of the matlab function to be optimized in the form
%               [objective, constraints] = func(x,consts)
%  x_init  :   the vector of initial parameter values                  (n x 1)
%  x_lb    :   lower bound on permissible values of the variables, x   (n x 1)
%  x_ub    :   upper bound on permissible values of the variables, x   (n x 1)
%  options :   options(1) = 1 means display intermediate results
%                           2 means display more intermediate results
%                           3, 4, ... even more and more intermediate results
%              options(2)  = tol_x  tolerance on convergence of design variables
%              options(3)  = tol_f  tolerance on convergence of objective
%              options(4)  = tol_g  tolerance on convergence of constraints
%              options(5)  = max_evals limit on number of function evaluations
%              options(6)  = penalty  on constraint violations
%              options(7)  = exponent on constraint violations
%              options(8)  = max number of function eval's in est. of mean f(x)
%              options(9)  = desired accuracy of mean f (as a c.o.v.)
%              options(10)  = 1  means stop when solution is feasible   
%  consts  :   an optional vector of constats to be passed to func(x,consts)
%
% OUTPUT
% ======
%  x_opt   :   a set of design variables at or near the optimal value
%  f_opt   :   the objective associated with the near-optimal design variables
%  g_opt   :   the constraints associated with the near-optimal design variables
% cvg_hst  :   record of x, f, g, function_count, and convergence criteria

 
% Sheela Belur V, An Optimized Step Size Random Search ,
% Computer Methods in Applied Mechanics & Eng'g, Vol 19, 99-106, 1979
%
% S.Rao, Optimization Theory and Applications, 2nd ed, John Wiley, 1984
%
% modifications by:
% H.P. Gavin , Civil & Env'ntl Eng'g, Duke Univ.

% ``Consider everything. Keep the good. Avoid evil whenever you notice it.''

BX = 1;

% handle missing arguments
 if nargin<3                      % default values
   x_lb = -1.0e2*abs(x_init);
   x_ub =  1.0e2*abs(x_init);
 end
 if nargin<5, options = []; end   % use default option values
 options = optim_options(options);
 if nargin<6, consts = 1.0; end;

 msglev    = options(1);         % 1: display intermediate values
 tol_x     = options(2);         % parameter  tolerance
 tol_f     = options(3);         % objective fctn  tolerance
 tol_g     = options(4);         % constraint tolerance
 max_evals = options(5);         % maximum number of function evaluations
 penalty   = options(6);         % penalty  on constraint violations
 q         = options(7);         % penalty  function exponent
 find_feas = options(10);        % stop once a feasible solution is found

 x_lb=x_lb(:); x_ub=x_ub(:); x_init=x_init(:);	% make column vectors

 n = length(x_init);                    % number of variables

 function_count = 0;			% number of function evaluations
 iteration = 1;
 cvg_f = 1;
 cvg_hst = NaN(n+5,max_evals);	        % convergence history
 fa = zeros(4,1);			% augmented cost function

 s0 = (x_lb+x_ub)./(x_lb-x_ub);  s1 = 2./(x_ub-x_lb);  % scale to [-1:+1]
 x_init = s0+s1.*x_init;                % scale x to [-1:+1]
 x_init = min(max(x_init,-1.0),1.0);

 i3 = 1e-3*eye(3);               % for regularization 

% Evaluate the initial guess and the range of allowable parameter variation
 if any(x_ub <= x_lb)
   disp('error: x_ub can not less than or equal to x_lb for any parameter');
   x_opt = (x_init-s0)./s1;
   f_opt = (sqrt(5)-1)/2;
   g_opt = 1;
   return
 else

   [fv,gv,x_init,cJ,nAvg] = avg_cov_func(func,x_init,s0,s1,options, consts,1);
   g_max = max(gv);
   function_count = function_count + nAvg;

   m = length(gv);                 % number of constraints

   if prod(size(fv)) ~= 1
    error('the objective computed by your analysis function must be a scalar')
   end
   if size(gv,1) == 1 && size(gv,2) > 1
    error('the constraints computped by your analysis function must be a column vector')
   end

 end

 if ( msglev > 2 )
   [f_min,f_max]=plot_surface(func,(x_init-s0)./s1,x_lb,x_ub,options,consts,103);
 end

% the standard deviation of fractional random perturbations  between 0 and 1
% larger  sigma ... more "broad exploring"
% smaller sigma ... more "narrow focusing"
 alpha = 1.00;   % ??? 0.72 ??? 
 sigma = 0.20;
 phi   = 0.85;

 if msglev, tic; end

% analyze initial guess
 x1 = x_init;
 [fa(1),g1,x1,c1,nAvg] = avg_cov_func(func,x1,s0,s1,options,consts,BX);
 function_count = function_count + nAvg;
 
% initialize optimal values 
 f_opt = fa(1);  x_opt = x1;  g_opt = g1;

% save the initial guess
 cvg_hst(:,iteration)=[(x_opt-s0)./s1;f_opt;max(g_opt);function_count;sigma;1];

 if msglev, more off; end

 x4 = x1; g4 = g1; fa(4) = fa(1);

 last_update = function_count;

%Gradient = randn(n,1);                  % random initial search direction

 while ( function_count < max_evals )		% === main loop ===

   r = sigma*randn(n,1);			% random search direction
%  r = -0.1*sigma * ( Gradient/norm(Gradient) + 0.0*rand(n,1) ); % SPSA search dir.

   a2 = box_constraint(x1,r,n);
   x2 = x1 + alpha*a2*r;			% +1 random perturbation

   [fa(2),g2,x2,c2,nAvg] = avg_cov_func(func,x2,s0,s1,options,consts,BX);
   function_count = function_count + nAvg;

%  % an SPSA-type gradient estimation for next search direction ????
%  Gradient = ( fa(2) - fa(1) ) ./ r ;

   if ( fa(2) < fa(1) ),  step = +2;  else  step = -1;  end
   a3 = box_constraint(x1,step*r,n);
   x3 = x1 + a3*step*r;				% +2 or -1 random purturbation
   [fa(3),g3,x3,c3,nAvg] = avg_cov_func(func,x3,s0,s1,options,consts,BX);
   function_count = function_count + nAvg;

   abc = [[0,0,1; 0.5*(alpha*a2)^2,alpha*a2,1; 0.5*(a3*step)^2,a3*step,1] + i3] \ fa(1:3);

   a = abc(1);   b = abc(2);   c = abc(3); 

   quad_update = 0;
   if a > 0 				% curvature is positive!
      d = -b/a;                         % try to go to zero-slope point
      a4 = box_constraint(x1,d*r,n);
      x4 = x1 + a4*d*r;
      [fa(4),g4,x4,c4,nAvg] = avg_cov_func(func,x4,s0,s1,options,consts,BX);
      function_count = function_count + nAvg;
   end
   if ( msglev > 2 )		   % plot values on the surface ----
     x1 = (x1-s0)./s1; x2 = (x2-s0)./s1; x3 = (x3-s0)./s1; x4 = (x4-s0)./s1;
     f_offset = (f_max-f_min)/100;
     figure(103)
       ii = options(11);  jj = options(12);
       plot3([x2(ii),x3(ii)],[x2(jj),x3(jj)],fa([2,3]) + f_offset,  '-or', 'MarkerSize', 6 );
       if quad_update
         plot3([x2(ii),x4(ii)],[x2(jj),x4(jj)],fa([2,4]) + f_offset, '-or', 'MarkerSize', 6 );
       end
%      plot3([x1(ii),x2(ii),x3(ii),x4(ii)],[x1(jj),x2(jj),x3(jj),x4(jj)], ...
%                                       fa + f_offset, 'or', 'MarkerSize', 4); 
       drawnow
     x1 = s0+s1.*x1; x2 = s0+s1.*x2; x3 = s0+s1.*x3; x4 = s0+s1.*x4; %pause(0.5)
   end                          % ---- plot values on the surface

   if ( msglev > 3 && a > 0 )
     figure(104);
       clf
       xx = [min([0,a3*step,a4*d])-0.5:0.1:max([alpha*a2,a3*step,a4*d])+0.5];
       ff = 0.5*a*xx.^2 + b*xx + c;
       hold on
       plot(sigma*xx,ff,'-', 'LineWidth', 6 );
       plot(sigma*[0,alpha*a2,a3*step,a4*d], fa,'o', ...
                                      'MarkerSize',19, 'LineWidth', 6 )
       plot(sigma*a4*d,fa(4),'*', 'MarkerSize',19, 'LineWidth', 6 )
       hold off
       title(sprintf(' iter:  %d  |  eval:  %d  |  f_{opt} = %11.4e  |  \\sigma = %5.3f', ...
                                  iteration, function_count, min(fa), sigma ));
       xlabel('change in paramter values, \Delta x')
       ylabel('augmented objective function, f_{aug}(x)')
       grid on;  axis('tight');  drawnow;   
   end
   if ( msglev > 3 && a < 0 )
      figure(104);
        clf
        plot( sigma*[0,alpha*a2,a3*step], fa(1:3),'o' , ...
                                      'MarkerSize',19, 'LineWidth',4 )
          title(sprintf(' iter:  %d  |  eval:  %d  |  f_{opt} = %11.4e  |  \\sigma = %5.3f', ...
                                  iteration, function_count, min(fa), sigma ));
          grid on;  axis('tight');  drawnow;    
   end
   
   [fa(1),i_min]=min(fa);               % find the lowest fa of 4 evaluations

   if ( i_min == 2 ), 	x1=x2; g1=g2; c1=c2; end % and the corresponding x,g,c
   if ( i_min == 3 ),	x1=x3; g1=g3; c1=c3; end
   if ( i_min == 4 ),	x1=x4; g1=g4; c1=c4; quad_update = 1; end

   if ( i_min > 1 )  % slightly reduce the scope of the search and continue
%     sigma = sigma / 1.1;
      sigma = sigma / function_count^(0.0200); 
%     sigma = sigma*(phi*exp(-function_count^2/0.2/max_evals^2) + (1-phi));
   end

   x1 = min(max(x1,-1),+1);        % ... keep x within limits [x_lb, x_ub] 
%

   if iteration > 1
      cvg_f = ( cvg_hst(n+1,iteration-1)/f_opt ) / ( function_count - cvg_hst(n+3,iteration-1));
   end

   if ( fa(1) < f_opt )  % update the optimal solution 

      x_opt = x1;   
      f_opt = fa(1);   
      g_opt = g1;
      last_update = function_count;
      iteration   = iteration + 1; % increment the iteration counter

      cvg_hst(:,iteration) = ...
            [ (x_opt-s0)./s1; f_opt; max(g_opt); function_count; sigma; cvg_f ];

      if msglev                     % display some results ----

         secs_left = round((max_evals-function_count)*toc/function_count);
         eta       = datestr(now+secs_left/3600/24,14);

         [max_g, idx_ub_g] = max(g_opt);
         home
         fprintf(1,' -+-+-+-+-+-+-+-+-+-+- ORSopt -+-+-+-+-+-+-+-+-+-+-+-+-+\n')
         fprintf(1,' iteration               = %5d', iteration );
         if max(g_opt) > tol_g 
            fprintf(1,'     !!! infeasible !!!\n');
         else
            fprintf(1,'     ***  feasible  ***\n');
         end
         fprintf(1,' function evaluations    = %5d  of  %5d  (%4.1f%%) \n', ...
                     function_count, max_evals, 100*function_count/max_evals );
         fprintf(1,' e.t.a.                  = %s\n', eta );
         fprintf(1,' objective               = %11.3e\n', f_opt );
         fprintf(1,' variables               = ');
         fprintf(1,'%11.3e', (x_opt-s0)./s1 );
         fprintf(1,'\n');
         fprintf(1,' max constraint          = %11.3e (%d) \n', max_g, idx_ub_g );
         fprintf(1,' Convergence Criterion F = %11.4e    tolF = %10.6f\n',cvg_f, tol_f);
         fprintf(1,' Convergence Criterion X = %11.4e    tolX = %10.6f\n',sigma, tol_x);
         fprintf(1,' c.o.v. of f_a  = %11.3e\n', c1 );
         fprintf(1,' -+-+-+-+-+-+-+-+-+-+- ORSopt -+-+-+-+-+-+-+-+-+-+-+-+-+\n')
         if quad_update
             fprintf(1,' line quadratic update successful\n');
             %%pause(1);
         end
      end                          % ---- display some results 


   end                             % ----- update the optimal point
%

   % check convergence criteria and display termination information -----

   if ( max(g_opt) < tol_g && find_feas ) % a feasible point has been found
      fprintf(1,' Woo Hoo!  Feasible solution found! ')
      fprintf(1,' ... and that is all we are asking for.\n')
      break;                       % ... and that's all we want
   end

   if ( iteration > 1 && ( (sigma<tol_x) || (cvg_f < tol_f) ) )
%        ||   (function_count-last_update>0.2*max_evals)))
      fprintf(1,' *** Woo-Hoo!  Converged solution found! \n')
      if ( sigma < tol_x )
        fprintf(1,' *** convergence in design variables \n');
      end
      if ( cvg_f < tol_f )
        fprintf(1,' *** convergence in objective \n');
      end
%     [f_opt,g_opt] = feval(func,(x_opt-s0)./s1,consts);
      if (max(g_opt) < tol_g )
        fprintf(1,' *** Woo-Hoo!  Converged solution is feasible! \n')
      else
        fprintf(1,' *** Boo-Hoo!  Converged solution is NOT feasible! \n')
        if max(g_opt) > tol_g
           fprintf(1,'     ... Increase options(6) and try, try again ...\n')
        else
           fprintf(1,'     ... Decrease options(6) and try, try again ...\n')
        end
      end
      if msglev, format, end
      break;
   end                         % ----- check convergence criteria

 end						% === main loop ===

 if function_count >= max_evals        % time-out --  
   if msglev
      disp([' Enough!!  Maximum number of function evaluations (', ...
              int2str(max_evals),') has been exceeded']);
      disp('     ... Increase tol_x (options(2)) or max_evals (options(5)) and try try again.')
      format
      %% current_options = options
   end
 end                                   % -- time-out

 x_init = (x_init-s0)./s1;              % scale back to original units
 x_opt  = (x_opt-s0)./s1;               % scale back to original units

 fprintf(1,' *** Completion  : %s (%s)\n', ...
 datestr(now,14), datestr(toc/3600/24,13));
 fprintf(1,' *** Objective   : %11.3e \n', f_opt);
 fprintf(1,' *** Variables   : \n');
 fprintf(1,'               x_init         x_lb     <     x_opt    <     x_ub\n')
 fprintf(1,'------------------------------------------------------------------\n');
 for ii = 1:nvars
   eqlb = ' ';
   equb = ' ';
   if ( x_opt(ii) < x_lb(ii) + tol_g + 10*eps )
     eqlb = '=';
   elseif ( x_opt(ii) > x_ub(ii) - tol_g - 10*eps )
     equb = '=';
   end
   fprintf(1,'x(%3u)  %12.5f   %12.5f %s %12.5f %s  %12.5f \n',ii,x_init(ii),x_lb(ii),eqlb,x_opt(ii),equb,x_ub(ii));
 end

 fprintf(1,' *** Constraints : \n');
 for jj=1:m
    binding = ' ';
    if ( g_opt(jj) > -tol_g )
      binding = ' ** binding ** ';
    end
    fprintf(1,'       g(%3u) = %12.5f  %s \n', jj, g_opt(jj), binding );
 end

 if msglev, format; end

 iteration = iteration + 1;            % save final iteration information 
 cvg_hst(:,iteration) = [x_opt; f_opt; max(g_opt); function_count; sigma;cvg_f];
 cvg_hst(n+5,[1,2]) = [1,1]*cvg_hst(n+5,3);
 cvg_hst = cvg_hst(:,1:iteration);

end
% ORSopt ======================================================================
% updated 2011-04-13, 2014-01-12, 2015-03-14, (pi day 03.14.15), 2015-03-26, 
% 2016-04-06, 2019-02-23, 2020-01-17, 2024-04-03
