function [ x , status ] =  GRS( f , varargin )

%function [ x , status ] = GRS( f , range , delta , MaxFeval )
%
% Apply the classical Golden Ratio Search for the minimization of the
% provided one-dimensional function f, which must have the following
% interface:
%
%   [ v , varargout ] = f( x )
%
% Input:
%
% - x is either a scalar real denoting the input of f(), or [] (empty).
%
% Output:
%
% - v (real, scalar): if x == [] this is the best known lower bound on
%   the global optimum of f() on the standard interval in which f() is
%   supposed to be minimised (see next). If x ~= [] then v = f(x).
%
% - g (real, either scalar or a [ 2 x 1 ] matrix denoting an interval) is
%   the first optional argument. This also depends on x. if x == [] then
%   g is a [ 2 x 1 ] matrix denoting the standard interval in which f()
%   is supposed to be minimised (into which v is the minimum). f() is
%   never called with x ~= [].
%
% The other [optional] input parameters are:
%
% - range: (either [ 2 x 1 ] real vector or [], default []): the range
%   in which the local minimum has to be seached; if range == [], the
%   default range point provided by f() is used.
%
% - delta (real scalar, default value 1e-6): the accuracy in the stopping
%   criterion: the algorithm is stopped when the diameter of the
%   restricted range is less than or equal to delta.
%
% - MaxFeval (integer scalar, default value 100): the maximum number of
%   function evaluations (hence, iterations will be not more than
%   MaxFeval - 2 because at each iteration one function evaluation is
%   performed, except in the first one when two are).
%
% Output:
%
% - x (real scalar): the best solution found so far.
%
% - status (string): a string describing the status of the algorithm at
%   termination
%
%   = 'optimal': the algorithm terminated having proven that x is a(n
%     approximately) optimal solution, i.e., the diameter of the
%     restricted range is less than or equal to delta.
%
%   = 'empty': the provided range is empty (x_- > x_+) and therefore
%     such is the optimization problem
%
%   = 'stopped': the algorithm terminated having exhausted the maximum
%     number of iterations: x is the best solution found so far, but not
%     necessarily the optimal one
%
%   = 'error': the algorithm found a numerical error that prevents it from
%     continuing optimization
%
%{
 =======================================
 Author: Antonio Frangioni
 Date: 09-08-21
 Version 0.10
 Copyright Antonio Frangioni
 =======================================
%}

Plotg = 1;
% 0 = nothing is plotted
% 1 = the function value / gap are plotted
% 2 = the function and the test points used are plotted

Interactive = true;  % if we pause at every iteration

% reading and checking input- - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if ~ isa( f , 'function_handle' )
   error( 'f not a function' );
end

if isempty( varargin ) || isempty( varargin{ 1 } )
   [ fStar , range ] = f( [] );
else
   fStar = - Inf;  % if the range is not the standard one, we can't trust
		   % the standard global minima
   range = varargin{ 1 };
end

if ~ isreal( range )
   error( 'range not a real vector' );
end

if ( size( range , 1 ) ~= 1 ) || ( size( range , 2 ) ~= 2 )
   error( 'range is not a [ 1 x 2 ] vector' );
end

xm = range( 1 );  % x_-
xp = range( 2 );  % x_+

if xm > xp
   fprintf( 'range is empty\n' );
   x = NaN;
   status = 'empty';
   return;
end

if length( varargin ) > 1
   delta = varargin{ 2 };
   if ~ isreal( delta ) || ~ isscalar( delta )
      error( 'delta is not a real scalar' );
   end
else
   delta = 1e-6;
end

if length( varargin ) > 2
   MaxFeval = round( varargin{ 3 } );
   if ~ isscalar( MaxFeval )
      error( 'MaxFeval is not an integer scalar' );
   end
   if MaxFeval < 2
      error( 'at least two function computations are required' );
   end
else
   MaxFeval = 100;
end

% initializations - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

r = ( sqrt( 5 ) - 1 ) / 2;

xmp = xm + ( 1 - r ) * ( xp - xm );  % x'_-
xpp = xm + r * ( xp - xm );          % x'_+

fxmp = f( xmp );                     % f( x'_- )
fxpp = f( xpp );                     % f( x'_+ )

feval = 2;

if fxmp <= fxpp
   fx = fxmp;
else
   fx = fxpp;
end

status = 'optimal';

if Plotg
   gap = [];
end

fprintf( 'Golden ratio search\n');
if fStar > - Inf
   fprintf( 'feval\trel gap\t\tx_-\t\tx_+\n');
else
   fprintf( 'feval\tfbest\t\tx_-\t\tx_+\n');
end

% main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

while xp - xm > delta

    % output statistics - - - - - - - - - - - - - - - - - - - - - - - - - -

    if fStar > - Inf
       gapk = ( fx - fStar ) / max( [ abs( fStar ) 1 ] );

       if Plotg == 1
	  gap( end + 1 ) = gapk;
	  semilogy( gap , 'Color' , 'k' , 'LineWidth' , 2 );
	  xlim( [ 0 35 ] );
	  ylim( [ 1e-15 inf ] );
	  ax = gca;
	  ax.FontSize = 16;
	  ax.Position = [ 0.03 0.07 0.95 0.92 ];
	  ax.Toolbar.Visible = 'off';
       end
    else
       gapk = fx;
    end
    fprintf( '%4d\t%1.4e\t%1.8e\t%1.8e\n' , feval , gapk , xm , xp );

    if Plotg == 2
       xbot = xm - ( xp - xm ) / 20;
       xtop = xp + ( xp - xm ) / 20;

       warning( 'off' , 'all' );
       fplot( @(x) f( x ) , [ xbot xtop ] , 'Color' , 'k' , ...
	      'LineWidth' , 1 );

       xlim( [ xbot xtop ] );
       yticks( [] );
       ax = gca;
       ax.FontSize = 16;
       ax.Position = [ 0.025 0.05 0.95 0.95 ];
       ax.Toolbar.Visible = 'off';
       xticks( [ xbot xm xmp xpp xp xtop ] );
       xticklabels( { num2str( xbot , '%1.1g' ) , 'x-' , 'x''-' , ...
		    'x''+' , 'x+' , num2str( xtop , '%1.1g' ) } );
       warning( 'on' , 'all' );
    end


    % stopping criteria - - - - - - - - - - - - - - - - - - - - - - - - - -

    if feval > MaxFeval
       status = 'stopped';
       break;
    end

    % main logic- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    if fxmp <= fxpp
       xp = xpp; xpp = xmp; xmp = xm + ( 1 - r ) * ( xp - xm );
       fxpp = fxmp; fx = fxmp; fxmp = f( xmp );
    else
       xm = xmp; xmp = xpp; xpp = xm + r * ( xp - xm );
       fxmp = fxpp; fx = fxpp; fxpp = f( xpp );
    end

    feval = feval + 1;

    if Interactive
       pause;
    end

    % iterate - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

end

% end of main loop- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% select final answer

if fxmp <= fxpp
   x = xmp;
else
   x = xpp;
end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

end  % the end- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -