% ch13_golden_section.m
% One-dimensional minimization by golden-section search.
% Based on the Chapter 13 example C13_F01B.m.
%
% Run this script. Change the objective and settings below for a new problem.
% The objective must return a finite real scalar and be unimodal on [a,b]:
% it decreases toward one minimum and increases beyond it. Derivatives are
% not needed. This method does not establish global optimality when the
% interval contains several local minima.

clear; close all; clc;

%% Problem and settings
f = @(x) (x - 2).^2 + 1;  % Example: minimum at x = 2, f = 1
a = -1.0;                 % Left endpoint
b =  4.0;                 % Right endpoint
tol_x = 1e-6;             % Absolute tolerance on final interval width
max_iter = 200;           % Maximum number of interval reductions
make_plot = true;         % Set false to suppress the function plot

%% Minimize
[xmin, fmin, info] = golden_section(f, a, b, tol_x, max_iter);

%% Results
fprintf('\nGolden-section search\n');
fprintf('Estimated minimum: x = %.12g\n', xmin);
fprintf('Objective value:   f = %.12g\n', fmin);
fprintf('Final interval:    [%.12g, %.12g]\n', info.a, info.b);
fprintf('Interval width:    %.6e\n', info.b - info.a);
fprintf('Iterations:        %d\n', info.iterations);
fprintf('Function calls:    %d\n', info.evaluations);
fprintf('Status:            %s\n', info.status);

%% Optional plot over the original search interval
if make_plot
    x_plot = linspace(a, b, 400);
    % Evaluate one point at a time: f need not accept a vector input.
    y_plot = arrayfun(f, x_plot);

    figure('Color', 'w');
    plot(x_plot, y_plot, 'k-', 'LineWidth', 1.5);
    hold on;
    plot(xmin, fmin, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 7);
    xlabel('x', 'FontSize', 14);
    ylabel('f(x)', 'FontSize', 14);
    title('Golden-section search', 'FontSize', 14, 'FontWeight', 'normal');
    legend('f(x)', 'Estimated minimum', 'Location', 'best');
    set(gca, 'FontSize', 14);
    xlim([a b]);
    grid on; box on;
    % Plot evaluations are separate from the search function-call count.
end

%% Reusable numerical routine
function [xmin, fmin, info] = golden_section(f, a, b, tol_x, max_iter)
% GOLDEN_SECTION Minimize a unimodal scalar function on a finite interval.
% Returns the final interval midpoint and its function value.
% info.converged means that the interval-width tolerance was reached.
% The interval is a numerical localization, not a rigorous error bound:
% rounded function values can prevent reliable comparisons near a minimum.
% Choose tol_x to suit the physical scale and objective's numerical accuracy.
%
% To use this routine from other programs, copy it and checked_value into
% a separate file named golden_section.m.

    validateattributes(a, {'double'}, {'real','scalar','finite'});
    validateattributes(b, {'double'}, {'real','scalar','finite','>',a});
    validateattributes(tol_x, {'double'}, ...
        {'real','scalar','finite','positive'});
    validateattributes(max_iter, {'double'}, ...
        {'real','scalar','finite','integer','nonnegative'});
    if ~isa(f, 'function_handle') || ~isfinite(b - a)
        error('Supply a function handle and an interval of finite width.');
    end

    rho = (sqrt(5) - 1) / 2;
    x1 = b - rho*(b - a);
    x2 = a + rho*(b - a);
    f1 = checked_value(f, x1);
    f2 = checked_value(f, x2);
    neval = 2;
    iter = 0;
    status = 'Maximum iteration count reached.';

    while (b - a > tol_x) && (iter < max_iter)
        if ~(a < x1 && x1 < x2 && x2 < b)
            status = 'Stopped: interval points cannot be resolved in floating point.';
            break;
        end

        % Retain the better side and reuse one point and its function value.
        % Each interval reduction requires only one new function evaluation.
        if f1 < f2
            b = x2;
            x2 = x1;
            f2 = f1;
            x1 = b - rho*(b - a);
            f1 = checked_value(f, x1);
        else
            a = x1;
            x1 = x2;
            f1 = f2;
            x2 = a + rho*(b - a);
            f2 = checked_value(f, x2);
        end
        iter = iter + 1;
        neval = neval + 1;
    end

    converged = (b - a <= tol_x);
    if converged
        status = 'Interval-width tolerance reached.';
    end

    xmin = a + (b - a)/2;
    fmin = checked_value(f, xmin);
    info = struct('a', a, 'b', b, 'iterations', iter, ...
        'evaluations', neval + 1, 'converged', converged, 'status', status);
end

function value = checked_value(f, x)
% Reject invalid objective values rather than returning a misleading result.
    value = f(x);
    if ~isnumeric(value) || ~isscalar(value) || ~isreal(value) || ~isfinite(value)
        error('The objective must return a finite real scalar at x = %.16g.', x);
    end
end