function [A B]=genalg(error, population_size, mutation, no_of_generations)

%% Genetic Algorithms function
% This function does the selection, mating, crossover, mutation and
% insertion (the fitness is evaluated at another function called fitness)
% [A B] = genalg()
% A a vector containing evaluated 'a' parameters
% B a vector containing evaluated 'b' parameters
%
% Author: Wesam ELSHAMY
% wesamelshamy@ieee.org
% Electrical Engineering Dept., Cairo University, Egypt

global Na Nb 
errortolerance = error;
popsize = population_size;
mutationrate = mutation;
count = 0;

%% Creating initial population
pop = 1*rand(popsize, Na + Nb + 1);
for x = 1 : popsize
    pop(x, Na + Nb + 1) = fitness(pop(x, :));
end
pop = sortrows(pop, Na + Nb + 1);                       % Sorting initial population

%% Generations loop
while pop(1, Na + Nb + 1) > errortolerance && count < no_of_generations
    gene = zeros(popsize - 1, Na + Nb + 1);
    
    %% ----- Crossover -----
    for x = 1 : popsize - 2
        new1 = pop(x + 1, :) + 0.6*rand*(pop(x + 2) - pop(x + 1,:));
        new2 = pop(x + 2, :) + 0.6*rand*(pop(x + 1) - pop(x + 2,:));
        crosspoint = round((Na + Nb - 3)*rand + 2);
        gene(x, :) = [new1(1 : crosspoint) new2(crosspoint + 1 : Na + Nb) 0];
        gene(x + 1, :) = [new2(1 : crosspoint) new1(crosspoint + 1 : Na + Nb) 0];
    end
        
    %% ----- Mutation -----
    no_of_mutations = round(mutationrate*popsize*(Na + Nb));
    mutation_column = round((Na + Nb - 1)*(rand(1, no_of_mutations)) + 1);
    mutation_row = round((popsize - 3)*(rand(1, no_of_mutations)) + 2);    
    for x = 1 : no_of_mutations
        gene(mutation_row(x), mutation_column(x)) = 2*rand*gene(mutation_row(x), mutation_column(x));         % Mutating a random bit
    end
    
    %% Fitness evaluation
    for x = 1 : popsize - 1
        gene(x, Na + Nb + 1) = fitness(gene(x, :));
    end
    
    %% Reeinserting the new citizens          
    temppop = sortrows([gene; pop], Na + Nb + 1);
    pop = temppop(1 : popsize, :);
    
    count = count + 1;
end

% Returned parameters

a = pop(1, 1 : Na);
b = pop(1, Na + 1 : Na + Nb);
A = zeros(1, size(a, 2));
B = zeros(1, size(b, 2));
for x = 1 : size(a, 2)
    A(x) = a(size(a, 2) - x + 1);
end
for x = 1 : size(b, 2)
    B(x) = b(size(b, 2) - x + 1);
end