%% Script File:  Optimization and GA
% Changes Hills.m to animate the results in a contour diagram...

clear

% Change for graphics:
[X,Y]=meshgrid(0:0.1:10);
Z=X.*sin(4*X)+1.1*Y.*sin(2*Y);
%[Con,h]=contour(X,Y,Z);

% Stopping criteria
maxit=50;  %Max number of iterations
mincost=-99999999;  %Minimum cost

% Parameters:
popsize=8;

% Initialize the population:
Pop=10*rand(popsize,2);
figure(1)
[Con,h]=contour(X,Y,Z);
hold on
plot(Pop(:,1),Pop(:,2),'k*');
hold off
MM(1)=getframe;


mutrate=0.001;  %Mutation rate
popKept=0.5;  %Fraction of the population to keep
keep=floor(popKept*popsize);  %How many individuals are kept
M=ceil((popsize-keep)/2); % number of matings; 2 mates create 2 offspring
crossprob=round(rand(maxit,1)); %0= x-coord, 1= y-coord
nmut=ceil((popsize-1)*2*mutrate);  %Number of mutations 

% Fitness function:
ff='testfunction';  %The fitness function is in the file testfunction.m

% Probability distribution (won't change in this example)
probs=(keep:-1:1)/sum(1:keep);  %Probability is rank ordering

%% Initialize the population:

cost=feval(ff,Pop);  %Initial costs
[cost,idx]=sort(cost);  % Default sort is from small to large
Pop=Pop(idx,:);  
minc(1)=min(cost);  %Minimum cost, for plotting later
meanc(1)=mean(cost);  %Mean cost for this population (for plotting later)


%% Main loop
iga=0;
while iga<maxit
    
    iga=iga+1;
    
    % Pair up and mate:
    ma=RandChooseN(probs,M);
    pa=RandChooseN(probs,M);

    % Set up crossover and mutation:
    idx2=keep+1:popsize;
    beta=rand;
    PopMa=Pop(ma,:);  PopPa=Pop(pa,:);
    if crossprob(iga)==0  %Crossover the x-coordinate
        for j=1:M
            Pop(idx2(2*j-1),1)=(1-beta)*PopMa(j,1)+beta*PopPa(j,1);
            Pop(idx2(2*j-1),2)=PopMa(j,2);
            Pop(idx2(2*j),1)=(1-beta)*PopPa(j,1)+beta*PopMa(j,1);
            Pop(idx2(2*j),2)=PopPa(j,2);
        end
    else  %Crossover the y-coordinate
        for j=1:M
            Pop(idx2(2*j-1),1)=PopMa(j,1);
            Pop(idx2(2*j-1),2)=(1-beta)*PopMa(j,2)+beta*PopPa(j,2);
            Pop(idx2(2*j),1)=PopPa(j,1);
            Pop(idx2(2*j),2)=(1-beta)*PopPa(j,2)+beta*PopMa(j,2);
        end
    end

    % Mutation

    mrow=sort(ceil(rand(1,nmut)*(popsize-1))+1);
    mcol=ceil(rand(1,nmut)*2);
    for ii=1:nmut
        Pop(mrow(ii),mcol(ii))=10*rand;
    end 
    
    % New cost, set up for next iteration:
    cost=feval(ff,Pop);  %Initial costs
    [cost,idx]=sort(cost);  % Default sort is from small to large
    Pop=Pop(idx,:);  
    minc(iga+1)=min(cost);  %Minimum cost, for plotting later
    meanc(iga+1)=mean(cost);  %Mean cost for this population (for plotting later)


    % Stopping criteria
    if iga>maxit || cost(1)<mincost
        break
    end
    
    % GRAPHICS:
    figure(1)
    [Con,h]=contour(X,Y,Z);
    hold on
    plot(Pop(:,1),Pop(:,2),'k*');
    if mod(iga,10)==0
        hold off
    end
    MM(iga+1)=getframe;
    
    
end % End of the while loop

%% Display the results
figure(2)
iters=0:length(minc)-1;
plot(iters,minc,iters,meanc,'-');
xlabel('generation');ylabel('cost');