% ch13_genetic_algorithm.m
% Companion code for Chapter 13: Optimization and Data Fitting.
%
% Bounded minimization by a real-valued genetic algorithm.
% Change the objective, bounds, and search settings below. Each objective
% call receives a column vector. No derivatives are required.
% Tournament selection compares objective values directly (smaller is better).
% Repeat with other seeds; global optimality is not guaranteed.

clear; close all; clc;

%% Problem and settings
f = @(p) sin(3*p(1))*cos(3*p(2)) + (p(1)-1)^2 + (p(2)+2)^2;
lower = [-3; -3];
upper = [ 3;  3];
population_size = 100;
max_generations = 200;
crossover_probability = 0.90;
mutation_probability = 0.30;  % Probability of mutating each child
mutation_sigma = 0.04*(upper-lower); % Gaussian displacement per coordinate
tournament_size = 3;
seed = 42;
make_plot = true;
plot_x = [-1, 3];
plot_y = [-3, 0];

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

%% Minimize
[pbest,fbest,info] = genetic_search(f,lower,upper,population_size, ...
    max_generations,crossover_probability,mutation_probability, ...
    mutation_sigma,tournament_size);
fprintf('\nGenetic algorithm\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),200), ...
                   linspace(plot_y(1),plot_y(2),200));
    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('Genetic algorithm','FontSize',14,'FontWeight','normal');
    legend([hc,hb],{'f(p_1,p_2)','Best point found'},'Location','best');
    set(gca,'FontSize',14); axis equal; xlim(plot_x); ylim(plot_y);
    grid on; box on;
    % Plot evaluations are excluded from the search function-call count.
end

%% Reusable numerical routine
function [pbest,fbest,info]=genetic_search(f,lower,upper,N,K,pc,pm,sigma,nt)
% One elite individual survives unchanged in each generation.
% Parents are selected by tournaments, sampling with replacement.
% Uniform crossover exchanges coordinates between two parents.
% Mutation applies a Gaussian displacement to all coordinates of a selected
% child. Bounds are enforced by clipping. Mutation widths remain fixed.
% Copy this routine and its helpers to genetic_search.m for reuse.
    validateattributes(lower,{'double'},{'real','vector','finite','nonempty'});
    validateattributes(upper,{'double'},{'real','vector','finite','nonempty'});
    validateattributes(sigma,{'double'},{'real','vector','finite','positive'});
    lower=lower(:); upper=upper(:); sigma=sigma(:); d=numel(lower);
    if numel(upper)~=d || numel(sigma)~=d || any(upper<=lower) || any(~isfinite(upper-lower))
        error('Supply matching bounds and mutation widths, with upper > lower.');
    end
    validateattributes(N,{'double'},{'scalar','integer','finite','>=',2});
    validateattributes(K,{'double'},{'scalar','integer','finite','positive'});
    validateattributes(nt,{'double'},{'scalar','integer','finite','positive'});
    validateattributes(pc,{'double'},{'scalar','real','finite','>=',0,'<=',1});
    validateattributes(pm,{'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
    neval=N;
    for k=1:K
        [~,idx]=min(values);
        next=zeros(d,N); next_values=zeros(1,N);
        next(:,1)=pop(:,idx); next_values(1)=values(idx);
        count=1;
        while count<N
            p1=select_parent(pop,values,nt);
            p2=select_parent(pop,values,nt);
            children=[p1,p2];
            if rand<pc
                mask=rand(d,1)<0.5;
                children(mask,1)=p2(mask);
                children(mask,2)=p1(mask);
            end
            for j=1:2
                if count==N
                    break;
                end
                child=children(:,j);
                if rand<pm
                    child=child+sigma.*randn(d,1);
                end
                child=max(min(child,upper),lower);
                count=count+1;
                next(:,count)=child;
                next_values(count)=checked_value(f,child);
                neval=neval+1;
            end
        end
        pop=next; values=next_values;
    end
    [fbest,idx]=min(values); pbest=pop(:,idx);
    info=struct('generations',K,'evaluations',neval);
end

function p=select_parent(pop,values,nt)
    candidates=randi(size(pop,2),nt,1);
    [~,j]=min(values(candidates));
    p=pop(:,candidates(j));
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