function [ a, b, it ] = opt_golden ( f_string, a, b, n, x_tol )

%*****************************************************************************80
%
%% opt_golden uses golden search to seek a minimizing point of f(x).
%
%  Discussion:
%
%    The user enters a formula for f(x), the starting interval [a,b],
%    the iteration limit n, and tolerances for the size of the final x interval.  
%
%    It is assumed that f(x) is unimodal in [a,b], that is, starting at
%    x=a, it is strictly decreasing until reaching some value x*, after
%    which it is strictly increasing until reaching x=b.
%
%    The program returns an estimate for the value of the minimizing x,
%    and the number of iterations.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%
%      [ x, it ] = opt_golden ( 'x^2+sin(4*x)', -1, 1, 25, 0.00001 )
%
%    or, if called with no arguments, it will request them:
%
%      [ x, it ] = opt_golden ( );
%      Enter function formula, like x^2: x^2+sin(4*x)
%      Enter left endpoint, a: -1
%      Enter right endpoint, b: 2
%      Enter number of iterations, n: 25
%      Enter x tolerance: 0.00001
%
%    The function is specified as a string which is either:
%    * a MATLAB expression using the argument 'x', 
%    * the name of an M-file followed by the argument '(x)'.
%
%    The string should not contain any spaces between symbols, except
%    when it is passed as a function argument in quotes.
%
%    It is not necessary to use the "dot" notation for expressions involving
%    '*', '/', or '^', but it doesn't hurt either.
%
%    Examples of function specifications:
%
%      x^2
%      x.^2
%      3/(x^4+5*x-6)
%      sin(7*x)*sqrt(x)/8
%      wiggle(x)     <-- where "wiggle.m" is a user-provided M file.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function to be minimized.
%
%    real a, b: the left and right endpoints.
%
%    integer n: the maximum number of iterations allowed.
%    n must be at least 1.
%
%    real x_tol: a tolerance for the width of the x interval.
%
%  Output:
%
%    real x: the estimate for the location of the minimizer of f(x).
%
%    integer it: the number of steps taken.
%

%
%  Get f, a, b, n, x_tol.
%
  if ( nargin < 1 )
    f_string = input ( 'Enter function formula, like x^2: ', 's' );
  end
  f_element = vector_to_element ( f_string );
  f = str2func ( [ '@(x)', '( ', f_element, ' )' ] );

  if ( nargin < 2 )
    a = input ( 'Enter left endpoint, a: ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end

  if ( nargin < 3 )
    b = input ( 'Enter right endpoint, b: ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end

  if ( nargin < 4 )
    n = input ( 'Enter iteration limit, n: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  if ( n < 1 )
    error ( 'N must be at least 1.' );
  end

  if ( nargin < 5 )
    x_tol = input ( 'Enter x tolerance, x_tol: ' );
  elseif ( ischar ( x_tol ) )
    x_tol = str2num ( x_tol );
  end

  if ( x_tol <= 0.0 )
    error ( 'x_tol must be positive.' );
  end

  [ a, b, it ] = golden_search ( f, a, b, n, x_tol );

  return
end
function [ a, b, it ] = golden_search ( f, a, b, n, x_tol )

%*****************************************************************************80
%
%% golden_search seeks a minimizer of a unimodal function in [a,b].
%
%  Discussion:
%
%    A unimodal function f(x) in [a,b] has a minimizer a <= c <= b such that
%    f(x) strictly decreases between a and c, and strictly increases between
%    c and b.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 August 2019
%
%  Author:
%
%    Original MATLAB code by J Nathan Kutz.
%    This version by John Burkardt.
%
%  Reference:
%
%    J Nathan Kutz,
%    Data Driven Modeling and Scientific Computation,
%    Oxford University Press, 2013,
%    ISBN: 978-0-19-966034-6.
%
%  Input:
%
%    function handle f: the function to be minimized.
%
%    real a, b: the left and right endpoints.
%
%    integer n: the maximum number of iterations allowed.
%    n must be at least 1.
%
%    real x_tol: a tolerance for the width of the x interval.
%
%  Output:
%
%    real a, b: the current interval containing the minimizer.
%
%    integer it: the number of steps taken.

  it = 0;
  g = ( - 1.0 + sqrt ( 5.0 ) ) / 2.0;

  x1 =         g   * a + ( 1.0 - g ) * b;
  x2 = ( 1.0 - g ) * a +         g   * b;
  f1 = f ( x1 );
  f2 = f ( x2 );

  for it = 1 : n

    if ( f1 < f2 )

      b = x2;
      x2 = x1;
      f2 = f1;
      x1 = g * a + ( 1.0 - g ) * b;
      f1 = f ( x1 );
    else

      a = x1;
      x1 = x2;
      f1 = f2;
      x2 = ( 1.0 - g ) * a + g * b;
      f2 = f ( x2 );
    end

    if ( abs ( b - a ) <= x_tol )
      return
    end

  end

  return
end
function f_element = vector_to_element ( f_vector )

%*****************************************************************************80
%
%% vector_to_element rewrites a MATLAB vector function elementwise.
%
%  Discussion:
%
%    In MATLAB, by default, the operators '*', '/', and '^' are interpreted
%    as linear algebraic operations.  In many cases, a user simply desired
%    the elementwise operation instead, which MATLAB indicates using a 
%    dotted notation: '.*', './', and '.^'.
%
%    This function converts a string to a form that uses dotted operators.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_vector, a string presumably representing a MATLAB expression.
%
%  Output:
%
%    string f_element, a copy of f_vector in which all operators have been
%    replaced by their dotted (elementwise) form.
%
  n_vector = length ( f_vector );
  f_element = [''];

  c_prev = '';
  for iv = 1 : n_vector
    c = f_vector(iv);
    if ( ( c == '*' || c == '/' || c == '^' ) && ( c_prev ~= '.' ) )
      f_element = [ f_element, ' .', c, ' ' ];
    else
      f_element = [ f_element, c ];
    end
    c_prev = c;
  end

  return
end

