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


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=[];

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
population(k,7) = immse(Y,Yind);
end


for iter = 1:max_gen
    %recombination
    temp = population;

    %crossover
    intermediateCrossover = intermediate_crossover(population,mu,populationSize);

    %mutation and evaluatuon mse
    offspring = mutation(temp,mu,populationSize,tau1,tau2,X,Y);

    %offspring selection
    population = select_population(offspring);

    iter=iter+1;
    errors=[errors population(1,7)]
    iter
    if population(1,7)<e
        break
    end

end

plot(errors)

function sorted = select_population(offspring)
    sorted = sortrows(offspring,7);
end

function off = mutation(pop, mu,populationSize,tau1,tau2,X,Y)
    x= 1;
    off = zeros(populationSize,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 = intermediate_crossover(population,mu,populationSize)
x=1;
offspring =[];
for i = 1:mu
    p1 = population(x,:);
    if x <= 1
        newRow = [p1(1)/2 p1(2)/2 p1(3)/2 p1(4)/2 p1(5)/2 p1(6)/2];
        temp = [offspring;newRow];
        offspring = temp;
    else
        p2 = population((x-1),:);
        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;

    end
    x = x+1;
    if x >= populationSize
        x=1;
    end

end
end