function [p_opt,f_opt,g_opt,cvg_hst]=ossrs(func,p_init,p_min,p_max,options)
% [p_opt,f_opt,g_opt,cvg_hst]=ossrs(func,p_init,p_min,p_max,options)
%
% OSSRS: Optimized Step Size Randomized Search
% Nonlienar optimization with inequality constraints using Random Search
%
%       minimizes f(p) such that g(p)<0 and p_min <= p_opt <= p_max.
%       f is a scalar cost function, p is a vector of parameters, and
%       g is a vector of constraints.
%
% INPUT
% ======
%  func    :   the name of the matlab function to be optimized in the form
%               [cost,constraints] = func(p)
%  p_init  :   the vector of initial parameter values ... a column vector
%  p_min   :   minimum permissible values of the parameters, p
%  p_max   :   maximum permissible values of the parameters, p
%  options :   options(1) = 1       means display intermediate results
%              options(2) = parameter  tolerence 
%              options(3) = cost fctn  tolerence 
%              options(4) = constraint tolerence 
%              options(5) = max_iter limit on function evaluations
%              options(6) = penalty  on constraint violations
%              options(7) = exponent on constraint violations
%              options(8) = 1  means stop when solution is feasible      
%
% OUTPUT
% ======
%  p_opt   :   a set of parameters at or near the optimal value
%  f_opt   :   the cost associated with the near-optimal parameters
%  g_opt   :   the constraints associated with the near-optimal parameters
% cvg_hst  :   convergence history
%
% 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.   1-7-06, 11-8-07, 6-8-08

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

% handle missing arguments
 if nargin<3                      % default values
   p_min = -1.0e6*abs(p_init);
   p_max =  1.0e6*abs(p_init);
 end
 if nargin<5, options = []; end   % use default option values
 options = optim_options(options);

 dsply     = options(1);         % 1: display intermediate values
 tol_p     = options(2);         % parameter  tolerance
 tol_f     = options(3);         % cost fctn  tolerance
 tol_g     = options(4);         % constraint tolerance
 max_iter  = options(5);         % maximum number of function evaluations
 penalty   = options(6);         % penalty  on constraint violations
 q         = options(7);
 find_feas = options(8);

 p_min=p_min(:); p_max=p_max(:); p_init=p_init(:);	% make column vectors
 n = length(p_init);
 p1 = p_init;

 function_count = 0;

 p1 = min(max(p_min,p1),p_max);			% apply constraints
 [J1,g1]=feval(func,p1);  			% initialize f at p
 f(1)  = J1 + penalty*sum(g1.*(g1>tol_g))^q;

 f_opt = f(1);
 p_opt = p1;

 function_count = function_count + 1;

 cvg_hst = zeros(n+3,max_iter);

% the standard deviation of fractional random perturbations  between 0 and 1
% larger  std_dev ... more "broad exploring"
% smaller std_dev ... more "narrow focusing"
 std_dev = 0.2;
 iteration = 0;

 if dsply  
    more off
 end

 if n == 2 & dsply == 2    % plot the surface	
  p4 = p1; f(4) = f(1);
  a1 = [p_min(1):(p_max(1)-p_min(1))/45:p_max(1)];
  a2 = [p_min(2):(p_max(1)-p_min(1))/50:p_max(2)];
  n1 = length(a1);
  n2 = length(a2);
  fex = zeros(n1,n2);
  for i=1:n1
      for j=1:n2
          a = [ a1(i) ; a2(j) ];
          [J,g] = feval(func,a);
          fex(i,j) = J + penalty*sum(g.*(g>tol_g))^q;
      end
  end
 end

 while ( function_count < max_iter )

   r = std_dev*(p_max-p_min).*randn(n,1);

   p2 = p1 - r;					% negative random perturbation
   p2 = min(max(p_min,p2),p_max);		% apply constraints
   [J2,g2]=feval(func,p2);  
   f(2) = J2 + penalty*sum(g2.*(g2>tol_g))^q;

   p3 = p1 + r;					% positive random purturbation
   p3 = min(max(p_min,p3),p_max);		% apply constraints
   [J3,g3]=feval(func,p3);  
   f(3) = J3 + penalty*sum(g3.*(g3>tol_g))^q;

   function_count = function_count + 2;

   quad_update = 0;
   a = f(2) - 2.0*f(1) + f(3);
   if a > 0 					% curvature is positive!
   
     lambda = -0.5*(f(3)-f(2))/a; 		% try to go to zero-slope point
     p4 = p1 + lambda*r;
     p4 = min(max(p_min,p4),p_max);		% apply constraints
     [J4,g4]=feval(func,p4);  
     f(4) = J4 + penalty*sum(g4.*(g4>tol_g))^q;

     function_count = function_count + 1;
   
   end

   if n == 2 & dsply == 2  % plot the surface
      figure(7)
      clf
      mesh(a1,a2,fex')
      hold on
      plot3([p1(1),p2(1),p3(1),p4(1)],[p1(2),p2(2),p3(2),p4(2)],f*1.05,'r*')
      hold off
       xlabel('p_1')
       ylabel('p_2')
       zlabel('f')
      drawnow
%     pause(1)
   end
   
   [f(1),i_min]=min(f);			% find the lowest f of 4 evaluations

   if ( i_min == 2 ), 	p1 = p2; g1 = g2; end	% and the corresponding p, g
   if ( i_min == 3 ),	p1 = p3; g1 = g3; end
   if ( i_min == 4 ),	p1 = p4; g1 = g4; quad_update = 1; end

   if f(1) < f_opt		 % the last 4 function eval's improved things

      iteration = iteration + 1; % increment the iteration counter


      f_opt = f(1);              % update the best solution found, so far.
      p_opt = p1;
      g_opt = g1;

      cvg_hst(1:n,iteration) = p_opt;
      cvg_hst(n+1,iteration) = f_opt;
      cvg_hst(n+2,iteration) = max(g_opt);
      cvg_hst(n+3,iteration) = std_dev;

      if dsply                   % display some results
         home
         clc
         fprintf(1,'-+-+-+-+-+-+-+-+-+-+- OSSRS -+-+-+-+-+-+-+-+-+-+-+-+-+-\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 \n', function_count );
         fprintf(1,'cost           = %11.3e\n', f_opt );
         fprintf(1,'parameters     = ');
         fprintf(1,'%11.3e', p_opt );
         fprintf(1,'\n');
         fprintf(1,'max constraint = %11.3e\n', max(g_opt) );
         fprintf(1,'stndrd deviatn = %11.3e\n', std_dev );
         fprintf(1,'-+-+-+-+-+-+-+-+-+-+- OSSRS -+-+-+-+-+-+-+-+-+-+-+-+-+-\n');
         if quad_update
            fprintf(1,'line quadratic update successful\n')
         end
      end

      if max(g_opt) < tol_g	% current minimum point is feasible enough
         if find_feas
            fprintf(1,' Woo Hoo!  Feasible solution found!\n')
            fprintf(1,' ... and that is all we are asking for.\n')
	    cvg_hst = cvg_hst(:,1:iteration);
            return
         end
         if std_dev < tol_p
            fprintf(1,' Woo Hoo!  Converged solution found!\n')
            fprintf(1,' ... and that is all we can ask for.\n')
	    cvg_hst = cvg_hst(:,1:iteration);
            return
         end
         std_dev = std_dev / 1.5;  % reduce the scope of the search
      end

   end
   cvg_hst = cvg_hst(:,1:iteration);
 end

% --------------------------------------------------------------------- OSSRS