% ch13_differential_evolution.m
% Companion code for Chapter 13: Optimization and Data Fitting.
%
% Bounded minimization by differential evolution: current-to-best/1/bin.
% Change the objective, bounds, and search settings below.
% Each objective call receives a column vector. No derivatives are needed.
% Repeat with different seeds; global optimality is not guaranteed.

clear; close all; clc;

%% Problem and settings
% Shifted Rastrigin function: exact minimum f = 0 at p(j) = 0.5.
f = @(p) sum((p-0.5).^2 + 10*(1-cos(2*pi*(p-0.5))));
lower = [-4.5; -4.5];
upper = [5.5; 5.5];
population_size = 50;
max_generations = 400;
F = 0.9;                  % Difference-vector scale
CR = 0.9;                 % Crossover probability per coordinate
seed = 3;
make_plot = true;
plot_x = [-1, 2];         % Display limits, independent of search bounds
plot_y = [-1, 2];         % Adjust for a different objective or solution

rng(seed, 'twister');
% Alternative: comment out the line above and uncomment the following.
% rng('shuffle');

%% Minimize
[pbest, fbest, info] = differential_evolution(f, lower, upper, ...
    population_size, max_generations, F, CR);
fprintf('\nDifferential evolution\n');
for j = 1:numel(pbest)
    fprintf('Best p(%d) = %.12g\n', j, pbest(j));
end
fprintf('Best objective value: %.12g\n', fbest);
fprintf('Generations:          %d\n', info.generations);
fprintf('Function calls:       %d\n', info.evaluations);
fprintf('Status: Search budget completed.\n');
fprintf('Repeat with other seeds; global optimality is not guaranteed.\n');

%% Optional contour plot
if make_plot && numel(pbest) == 2
    [X,Y] = meshgrid(linspace(plot_x(1),plot_x(2),250), ...
                     linspace(plot_y(1),plot_y(2),250));
    Z = arrayfun(@(x,y) f([x;y]), X, Y);
    figure('Color','w');
    [~,hc] = contour(X,Y,Z,25,'LineColor',[0.55 0.55 0.55]);
    hold on;
    hb = plot(pbest(1),pbest(2),'k+','MarkerSize',11,'LineWidth',1.5);
    xlabel('p_1','FontSize',14); ylabel('p_2','FontSize',14);
    title('Differential evolution','FontSize',14,'FontWeight','normal');
    legend([hc,hb],{'f(p_1,p_2)','Best point found'},'Location','eastoutside');
    set(gca,'FontSize',14); axis equal;
    xlim(plot_x); ylim(plot_y);
    grid on; box on;
    % Plot calls are excluded from the search evaluation count.
end

%% Reusable numerical routine
function [pbest,fbest,info] = differential_evolution(f,lower,upper,N,K,F,CR)
% Each generation uses a frozen population and its best member.
% Two distinct random donors exclude the target particle.
% At least one coordinate is taken from the mutant (binomial crossover).
% Out-of-bounds coordinates are resampled uniformly within their bounds.
% Greedy selection preserves each member unless its trial is no worse.
% Copy this routine and its helper to differential_evolution.m for reuse.
    validateattributes(lower,{'double'},{'real','vector','finite','nonempty'});
    validateattributes(upper,{'double'},{'real','vector','finite','nonempty'});
    lower=lower(:); upper=upper(:); d=numel(lower);
    if numel(upper)~=d || any(upper<=lower) || any(~isfinite(upper-lower))
        error('Supply matching finite bounds with upper > lower.');
    end
    validateattributes(N,{'double'},{'scalar','integer','finite','>=',4});
    validateattributes(K,{'double'},{'scalar','integer','finite','positive'});
    validateattributes(F,{'double'},{'scalar','real','finite','>',0,'<=',2});
    validateattributes(CR,{'double'},{'scalar','real','finite','>=',0,'<=',1});
    pop=repmat(lower,1,N)+repmat(upper-lower,1,N).*rand(d,N);
    values=zeros(1,N);
    for i=1:N
        values(i)=checked_value(f,pop(:,i));
    end
    for k=1:K
        [~,idx]=min(values); best=pop(:,idx);
        next=pop; next_values=values;
        for i=1:N
            candidates=[1:i-1 i+1:N];
            donors=candidates(randperm(N-1,2));
            mutant=pop(:,i)+F*(best-pop(:,i)) ...
                +F*(pop(:,donors(1))-pop(:,donors(2)));
            mask=rand(d,1)<CR;
            mask(randi(d))=true;
            trial=pop(:,i); trial(mask)=mutant(mask);
            outside=trial<lower | trial>upper | ~isfinite(trial);
            trial(outside)=lower(outside)+(upper(outside)-lower(outside)).*rand(sum(outside),1);
            value=checked_value(f,trial);
            if value<=values(i)
                next(:,i)=trial; next_values(i)=value;
            end
        end
        pop=next; values=next_values;
    end
    [fbest,idx]=min(values); pbest=pop(:,idx);
    info=struct('generations',K,'evaluations',N*(K+1));
end

function value=checked_value(f,p)
    value=f(p);
    if ~isnumeric(value) || ~isscalar(value) || ~isreal(value) || ~isfinite(value)
        error('The objective must return a finite real scalar.');
    end
end