function [ a, b, mx, my, k ] = optisection ( a, b, c, f, e, n )

%*****************************************************************************80
%
%% optisection() seeks a minimizer of a unimodal function in [a,b].
%
%  Discussion:
%
%    A strict unimodal function f(x) in [a,b] has a minimizer c 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:
%
%    02 May 2026
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Vladimir Kodnyanko,
%    Description and software implementation of optisection method
%    for minimizing unimodal functions.
%
%    Vladimir Kodnyanko,
%    Optimal section method for minimizing unimodal functions,
%    Computational Mathematics and Modeling,
%    Volume 36, pages 235-252, 16 October 2025.
%
%  Input:
%
%    real a, b: the left and right endpoints.
%
%    real c: the convex combination coefficient.  A suggested value is 0.2.
%
%    function handle f: the function to be minimized.
%
%    real e: a tolerance.
%
%    integer n: the maximum number of iterations allowed.
%    n must be at least 1.
%
%  Output:
%
%    real a, b: the current interval containing the minimizer.
%
%    real mx, my: the estimated minimizer, and the minimizing value.
%
%    integer k: the number of function evaluations.
%
  k = 0;
%
%  At the beginning of at iteration, a 
%  the minimizer is estimated to occur at a <= mx <= b, with value f(mx) = my.
%
  for it = 1 : n
%
%  Special handling of first step.
%
    if ( it == 1 )

      mx = ( a + b ) / 2.0;
      my = f ( mx );
      k = k + 1;
%
%  The current minimizer is at mx, with function value my.
%  A new test point px is determined, and its function value py compared to my.
%  Depending on the result, the value of mx may be updated.
%  The interval [a,b] is updated.
%
    else

      if ( b - mx <= mx - a )
        px = c * a + ( 1.0 - c ) * mx;
      else
        px = c * b + ( 1.0 - c ) * mx;
      end

      py = f ( px );
      k = k + 1;

      if ( py <= my )

        if ( px <= mx )
          b = mx;
        else
          a = mx;
        end

        mx = px;
        my = py;

      else

        if ( mx < px )
          b = px;
        else
          a = px;
        end

      end
%
%  Decide if we can stop.
%
      e0 = e * abs ( mx ) + 1.0e-10;
      if ( abs ( mx - 0.5 * ( a + b ) ) <= e0 - 0.5 * ( b - a )  )
        return
      end

    end

  end

  return
end

