maxGen = 50;
populationSize = 100;
mu = 100;
e = 0.2;
population = zeros(populationSize,7); %a b c aStd bStd cStd MSE


fileID = fopen('model26.txt');
C = textscan(fileID,'%f %f');
fclose(fileID);
asd = C(1);
asdd = C(2);
X = [asd{:}];
Y = [asdd{:}];

tau1 = 1 / sqrt(2 * populationSize);
tau2 = 1 / sqrt(2 * sqrt(populationSize));

errors=[];

%START Timer
tic

%Generation of the initial population
for k = 1:populationSize
population(k,1:3) = normrnd(0,1,[1,3]);
population(k,4:6) = abs(normrnd(0,1,[1,3]));
Yind = zeros(101,1);
for i=1:length(X)
    Yind(i)=func(population(k,1),population(k,2),population(k,3),X(i));
end
%MSE for each individual
population(k,7) = immse(Y,Yind);
end


for iter = 1:maxGen
    %recombination
    %tmpPop = population;
    tmpPop = intermediateCrossover(population,mu,populationSize);
    %tmpPop = discreteCrossover(population,mu,populationSize);

    %mutation and evaluatuon mse for each individual
    offspring = mutation(tmpPop,mu,populationSize,tau1,tau2,X,Y);

    %offspring selection
    %population = selectPopulation(offspring,populationSize);
    population = selectPopulation2(offspring,population,populationSize);
    iter=iter+1;
    errors=[errors population(1,7)]
    iter
    if population(1,7)<e
        break
    end

end

final = population(1,1:7);
disp("a: "+final(1)+", b: "+final(2)+", c: "+final(3)+", MSE: "+final(1,7))
YY = zeros(101,1);
for i=1:length(X)
    YY(i)=func(final(1,1),final(1,2),final(1,3),X(i));
end

t = toc
%
figure(1);
plot(X,Y)
title('Plot')
hold on
plot(X,YY)
legend('Real values','Computed')
hold off
figure(2);
title('MSE')
plot(errors)

%(mu,lambda)
function sorted = selectPopulation(offspring,populationSize)
    sorted = sortrows(offspring,7);
    sorted = sorted(1:populationSize,:);
end
%(mu+lambda)
function sorted = selectPopulation2(offspring,parents,populationSize)
    sorted = sortrows([parents ;offspring],7);
    sorted = sorted(1:populationSize,:);
end


function off = mutation(pop, mu,populationSize,tau1,tau2,X,Y)
    x= 1;
    off = zeros(mu,7);

    for i =1:mu
       p = pop(x,:);
       off(i,1:6) = [(p(1)+normrnd(0,p(4))) (p(2)+normrnd(0,p(5))) (p(3)+normrnd(0,p(6))) (p(4)*exp(normrnd(0,tau1))*exp(normrnd(0,tau2))) (p(5)*exp(normrnd(0,tau1))*exp(normrnd(0,tau2))) (p(6)*exp(normrnd(0,tau1))*exp(normrnd(0,tau2)))];
       x=x+1;
       if x>=populationSize
           x=1;
       end
    end
    for k = 1: length(pop)
        Yind = zeros(101,1);
        for i=1:length(X)
            Yind(i)=func(off(k,1),off(k,2),off(k,3),X(i));
        end
        off(k,7) = immse(Y,Yind);
    end

end

function f = func(a,b,c,x)
  f = a*(x.^2-b*cos(c*pi*x));
end

function offspring = intermediateCrossover(population,mu,populationSize)
    x=1;
    offspring =[];
    for i = 1:mu
        p1 = population(x,:);
        if x <= 1
            p2 = population((populationSize),:);
        else
            p2 = population((x-1),:);
        end
        newRow = [((p1(1)+p2(1))/2) ((p1(2)+p2(2))/2) ((p1(3)+p2(3))/2) ((p1(4)+p2(4))/2) ((p1(5)+p2(5))/2) ((p1(5)+p2(5))/2)];
        temp = [offspring;newRow];
        offspring = temp;
        x = x+1;
        if x >= populationSize
            x=1;
        end

    end
end

function offspring = discreteCrossover(population,mu,populationSize)
    x =1;
    offspring = [];
    for i=1:mu
         p1 = population(x,:);
         if x <= 1
            p2 = population(populationSize,:);
            randomSelect = [0 0 0 0 0 0];
            for g = 1:6
                if rand() > 0.5
                    randomSelect(g) = p1(g);
                else
                    randomSelect(g) = p2(g);
                end
            end
            temp = [offspring;randomSelect];
            offspring = temp;
        else
            p2 = population(x-1,:);
            randomSelect = [0 0 0 0 0 0];
            for g = 1:6
                if rand() > 0.5
                    randomSelect(g) = p1(g);
                else
                    randomSelect(g) = p2(g);
                end
            end
            temp = [offspring;randomSelect];
            offspring = temp;
         end
          x = x+1;
        if x >= populationSize
            x=1;
        end
        offspring;
    end
end