! Illustrate Golden Section line search to find a small interval
 around the minimum of a uni-modal/quasi-convex function;
! Keywords: Golden section, Line search, Fibonacci search;
procedure compf:
! A function for which we want to find the minimum;
  fof = (x-2)^2;
ENDprocedure

procedure goldensection:
! Golden section search to find the minimum of f(x)
   for L <= x  <= U,
 Inputs:
    L = lower limit on value of x,
    U = upper limit on value of x,
    gstol = tolerance for accuracy of x,
    fof = function to compute f(x),
;
   goldratio = (-1+5^(0.5))/2;
   grc = 1 - goldratio; ! Its complement;
   @WRITE(' Golden ratio= ', goldratio, @NEWLINE(1));
   a = L;
   b = U;
! Graphically, the function is evaluated at the two points c and d
  in the interval:  a---c--d---b.
  If f(c) < f(d), then we know the min is in [a, d],
                  else, the min is in [c, b].
   In either case, add a new point, and repeat.
   Using the Golden section/ratio to choose the location of the
   new point approximately minimizes the number of iterations
   needed to achieve a specified tolerance interval. The Golden
   section search has behavior almost as good as the more
   complicated Fibonacci search;
!   Choose initial points;
   c = a + (b - a) * grc;
   d = a + (b - a) * goldratio; 
   x = c;
   compf; ! Compute function value at c;
   fc = fof;
   x = d;
   compf; ! Compute function value at d;
   fd = fof;
   iter = 0;
!  Each pass reduces the interval of uncertainty by the "golden ratio", 0.618;
   @WHILE( @ABS(b-a) #gt# gstol:
      iter = iter + 1;
      @WRITE('  fc, fd = ', fc, ' ', fd, @NEWLINE(1));
      @IFC( fc #lt# fd:
            b = d; ! New interval is [a, d];
            d = c;
            fd = fc;
            c = a + grc * (b-a);
            x = c;
            compf;
            fc = fof;
        @else
            a = c; ! New interval is [c, b];
            c = d;
            fc = fd;
            d = b - grc * (b-a);
            x = d;
            compf;
            fd = fof;
           );
      @WRITE(' Interval(',iter,') = [',a,', ',b,']', @NEWLINE(1));
        );  
   gsxstar = (a+b)/2; ! Return an estimate;
ENDprocedure

CALC:
  @SET( 'TERSEO',1);    ! Output level (0:verb, 1:terse, 2:only errors, 3:none);
  @SET( 'OROUTE',1);    ! Route output immediately to the window line by line;
  L = 1; ! Lower bound on x;
  U = 6; ! Upper bound on x;
  gstol = 0.01;  ! Desired accuracy;
  goldensection;
  @WRITE(' The estimated xstar = ', gsxstar, @NEWLINE(1));
ENDCALC