% ch13_simulated_annealing.m
% Companion code for Chapter 13: Optimization and Data Fitting.
%
% One-dimensional minimization by simulated annealing on a finite interval.
% Change the objective, bounds, starting point, and search settings below.
% Uphill moves can help the search escape a local minimum. Cooling gradually
% reduces their acceptance probability. No derivatives are required.
% The output is the best point found, not proof of a global minimum.
% Repeat with different seeds and cooling schedules to assess reliability.

clear; close all; clc;

%% Problem and settings
f = @(x) x.^4 - 20*x.^2 + 24*x + 5;
a = -5.0;                 % Lower bound
b =  4.0;                 % Upper bound
x0 = 2.0;                 % Starting point

T0 = 100.0;                % Initial temperature (in objective-value units)
T_min = 1e-3;             % Stop cooling below this temperature
cooling = 0.95;           % Temperature multiplier, between 0 and 1
moves_per_temp = 100;     % Number of trial moves at each temperature
step_size = 0.5;          % Maximum absolute displacement of a trial move
max_levels = 1000;        % Safety limit on number of temperature levels
seed = 1;                % Change to repeat with a different random sequence
make_plot = true;

rng(seed, 'twister');

%% Minimize
[xbest, fbest, info] = simulated_annealing(f, a, b, x0, ...
    T0, T_min, cooling, moves_per_temp, step_size, max_levels);

%% Results
fprintf('\nSimulated annealing\n');
fprintf('Best point found:     x = %.12g\n', xbest);
fprintf('Best objective value: f = %.12g\n', fbest);
fprintf('Temperature levels:   %d\n', info.levels);
fprintf('Trial moves:          %d\n', info.trials);
fprintf('Function calls:       %d\n', info.evaluations);
fprintf('Accepted moves:       %.2f %%\n', 100*info.acceptance_fraction);
fprintf('Out-of-bounds trials: %d\n', info.outside);
fprintf('Temperature at stop:  %.6e\n', info.temperature);
fprintf('Status:               %s\n', info.status);
fprintf('Repeat with other seeds; global optimality is not guaranteed.\n');

%% Optional function plot
if make_plot
    xx = linspace(a, b, 500);
    yy = arrayfun(f, xx);
    figure('Color', 'w');
    plot(xx, yy, 'k-', 'LineWidth', 1.5);
    hold on;
    plot(xbest, fbest, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 7);
    xlabel('x', 'FontSize', 14);
    ylabel('f(x)', 'FontSize', 14);
    title('Simulated annealing', 'FontSize', 14, 'FontWeight', 'normal');
    legend('f(x)', 'Best point found', 'Location', 'best');
    set(gca, 'FontSize', 14);
    xlim([a b]); grid on; box on;
    % Plot evaluations are excluded from the reported function-call count.
end

%% Reusable numerical routine
function [xbest, fbest, info] = simulated_annealing(f, a, b, x0, ...
        T0, T_min, cooling, moves_per_temp, step_size, max_levels)
% Symmetric uniform proposals and Metropolis acceptance at each temperature.
% Out-of-bounds proposals are rejected, rather than clipped to a boundary.
% A rejected move leaves the current state unchanged.
% The best state is stored separately: the current state may move uphill.
% Copy this routine and its helper to simulated_annealing.m for reuse.

    validateattributes(a, {'double'}, {'real','scalar','finite'});
    validateattributes(b, {'double'}, {'real','scalar','finite','>',a});
    validateattributes(x0, {'double'}, {'real','scalar','finite','>=',a,'<=',b});
    validateattributes(T0, {'double'}, {'real','scalar','finite','positive'});
    validateattributes(T_min, {'double'}, ...
        {'real','scalar','finite','positive','<',T0});
    validateattributes(cooling, {'double'}, {'real','scalar','finite','>',0,'<',1});
    validateattributes(moves_per_temp, {'double'}, ...
        {'real','scalar','finite','integer','positive'});
    validateattributes(step_size, {'double'}, {'real','scalar','finite','positive'});
    validateattributes(max_levels, {'double'}, ...
        {'real','scalar','finite','integer','positive'});

    x = x0;
    fx = checked_value(f, x);
    xbest = x;
    fbest = fx;
    T = T0;
    levels = 0; trials = 0; accepted = 0; outside = 0; neval = 1;

    while T > T_min && levels < max_levels
        for j = 1:moves_per_temp
            trials = trials + 1;
            x_trial = x + step_size*(2*rand - 1);
            if ~isfinite(x_trial) || x_trial < a || x_trial > b
                outside = outside + 1;
                continue;
            end
            f_trial = checked_value(f, x_trial);
            neval = neval + 1;
            df = f_trial - fx;

            % Accept all downhill moves; accept uphill moves with exp(-df/T).
            % The logarithmic form avoids evaluating a very small exponential.
            if df <= 0 || log(rand) < -df/T
                x = x_trial;
                fx = f_trial;
                accepted = accepted + 1;
                if fx < fbest
                    xbest = x;
                    fbest = fx;
                end
            end
        end
        levels = levels + 1;
        T = cooling*T;
    end

    if T <= T_min
        status = 'Cooling schedule completed.';
    else
        status = 'Maximum temperature-level count reached.';
    end
    info = struct('levels', levels, 'trials', trials, ...
        'evaluations', neval, 'acceptance_fraction', accepted/trials, ...
        'outside', outside, 'temperature', T, 'status', status);
end

function value = checked_value(f, x)
    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