%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   A demo of the genetic algorithm using a 2D example  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   Prepared by Mr. Christopher Swierczewski            %%%
%%%   Feb 6, 2015 for AMATH301 U of Washington, Seattle   %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

set(0, 'defaultLineLineWidth', 2)
set(0, 'DefaultAxesFontSize', 16)

%% Animating the Genetic Algorithm
%

f = @(x,y) x.*cos(x).*sin(y);

% Plot the Objective Function
xl = 0; xr = 8;
yl = -5; yr = 5;
xpoints = xl:0.2:xr;
ypoints = yl:0.2:yr;
[X,Y] = meshgrid(xpoints,ypoints);
F = f(X,Y);

figure(1)
contour(X,Y,F,16)
title('Animation will begin in 5 seconds...')


%% Genetic Algorithm
%
% We want to find the (x,y)-points that minimize f
%
% (CHANGING "rng(1)" TO "rng(0)" WILL LEAD TO A DIFFERENT RESULT!)
%


rng(1);
num_generations = 4;
num_mutations_per_survivor = 5;
num_survivors_per_generation = 2;
total_mutations_per_generation = num_mutations_per_survivor * num_survivors_per_generation;



% initial (random) guesses
xavg = (xl+xr)/2;
yavg = (yl+yr)/2;
mutations_x = xavg + randn(num_mutations_per_survivor, 1);
mutations_y = yavg + randn(num_mutations_per_survivor, 1);
survivors_x = mutations_x;
survivors_y = mutations_y;

% (pause for a moment so we can full-screen the plot)
pause(5)

for generation_number = 1:num_generations
    % evaluate the objective function at each current survivors
    f_values = f(mutations_x, mutations_y);
    
    % sort the f_values from smallest to largest and grab the indices that
    % result in the smallest values of f
    [f_values_sorted, indices] = sort(f_values);
    
    % grab the num_survivors_per_generation best approximations
    best_indices = indices(1:num_survivors_per_generation);
    survivors_x = mutations_x(best_indices);
    survivors_y = mutations_y(best_indices);
    
    % generate the next set of mutations
    new_mutations_x = [];
    new_mutations_y = [];
    for n = 1:num_survivors_per_generation
        mutations_xn = survivors_x(n) + randn(num_mutations_per_survivor,1)/(generation_number);
        mutations_yn = survivors_y(n) + randn(num_mutations_per_survivor,1)/(generation_number);

        new_mutations_x = [new_mutations_x, mutations_xn];
        new_mutations_y = [new_mutations_y, mutations_yn];
    end
    
    %
    % plot the new mutations and survivors
    %
    contour(X,Y,F,16)
    hold on
    plot(mutations_x, mutations_y, 'k.', 'Markersize', 20)
    pause(1)
    plot(survivors_x, survivors_y, 'ro', 'Markersize', 10)
    pause(1)
    for m = 1:num_survivors_per_generation
        xm = survivors_x(m);
        ym = survivors_y(m);
        for n = 1:num_mutations_per_survivor
            x = new_mutations_x(n,m);
            y = new_mutations_y(n,m);
            x_line = [xm x];
            y_line = [ym y];
            plot(x_line, y_line, 'r-', 'Linewidth', 1)
            plot(x, y, 'r.', 'Markersize', 20)
        end
    end
    %
    % End plotting per step
    %
    pause(1)
    hold off
        
    % update
    mutations_x = reshape(new_mutations_x,[total_mutations_per_generation,1]);
    mutations_y = reshape(new_mutations_y,[total_mutations_per_generation,1]);
end

%% Plot final best guess
%

f_values = f(mutations_x, mutations_y);
[f_values_sorted, indices] = sort(f_values);
best_index = indices(1);
xbest = mutations_x(best_index);
ybest = mutations_y(best_index);
fbest = f_values_sorted(1);

hold off
contour(X,Y,F,16)
hold on
plot(xbest, ybest, 'r.', 'Markersize', 24);
title(['Minimum Value f = ' num2str(fbest)])


