function [best_fitness, best_MMRE, best_PRED_25, best_Sum_Ab_Res, best_Med_Ab_Res, best_SD_Ab_Res, best_number_of_removed_features] = ...
          binary_GA(population_size, elitism_population, number_of_generations, crossover_rate, mutation_rate, regression_method, X_train, Y_train, X_test, Y_test, dataset_name)

global report_progress_5; % report index of generation in GA
global report_progress_6; % report index of chromosome in GA

if strcmp(regression_method, 'SVR with RBF kernel')
    number_of_parameters = 3;
    range_parameters{1} = [80, 150];   %--> C
    range_parameters{2} = [0.01, 1]; %--> epsilon
    range_parameters{3} = [0.01, 1]; %--> gamma
elseif strcmp(regression_method, 'SVR with linear kernel')
    number_of_parameters = 2;
    range_parameters{1} = [80, 150];   %--> C
    range_parameters{2} = [0.01, 1]; %--> epsilon
elseif strcmp(regression_method, 'MLP')
    number_of_parameters = 4;
    range_parameters{1} = [3, 10];   %--> N
    range_parameters{2} = [10, 30]; % [100, 5000]; %--> E (epochs)
    range_parameters{3} = [0.01, 1]; %--> L
    range_parameters{4} = [0.01, 1]; %--> M
elseif strcmp(regression_method, 'M5P')
    number_of_parameters = 3;
    range_parameters{1} = [2, 20];   %--> I
    range_parameters{2} = [0, 1]; %--> P
    range_parameters{3} = [0, 1]; %--> S
end
number_of_features = size(X_train, 2);
number_of_bits_per_parameter = 20;
dimension = (number_of_parameters * number_of_bits_per_parameter) + number_of_features;

%%%%% initialize population:
chromosomes = zeros(population_size, dimension);
for chromosome_index = 1:population_size
    while true
        chromosomes(chromosome_index, :) = round(rand(1, dimension));
        selected_features = chromosomes(chromosome_index, end-number_of_features+1 : end);
        if ~sum(selected_features) == 0  %--> not all features are removed
            break
        end
    end
end

for generation = 1:number_of_generations
    if report_progress_5 == true; str = sprintf('generation #%d', generation); disp(str); end;
    fitness = zeros(population_size, 1); MMRE = zeros(population_size, 1);
    PRED_25 = zeros(population_size, 1); Sum_Ab_Res = zeros(population_size, 1);
    Med_Ab_Res = zeros(population_size, 1); SD_Ab_Res = zeros(population_size, 1);
    %%%%% process every chromosome:
    for chromosome_index = 1:population_size
        if report_progress_6 == true; str = sprintf('----- chromosome #%d', chromosome_index); disp(str); end;
        chromosome = chromosomes(chromosome_index, :);
        %%%%% extract parameters out of chromosome:
        parameters = zeros(number_of_parameters, 1);
        for parameter_index = 1:number_of_parameters
            d_binary = chromosome((parameter_index-1)*number_of_bits_per_parameter + 1 : ...
                                   parameter_index*number_of_bits_per_parameter);
            d_decimal = 0;
            for bit_index = 1:number_of_bits_per_parameter
                if d_binary(bit_index) == 1
                    d_decimal = d_decimal + (2^(bit_index-1));
                end
            end
            min_parameter = range_parameters{parameter_index}(1);
            max_parameter = range_parameters{parameter_index}(2);
            parameter = min_parameter + (((max_parameter - min_parameter) / (2^number_of_bits_per_parameter - 1)) * d_decimal);
            parameters(parameter_index) = parameter;
        end
        %%%%% extract selected features out of chromosome:
        selected_features = chromosome(end-number_of_features+1 : end);
        %%%%% evaluate the fitness of chromosomes:
        [fitness(chromosome_index), MMRE(chromosome_index), PRED_25(chromosome_index), ...
            Sum_Ab_Res(chromosome_index), Med_Ab_Res(chromosome_index), SD_Ab_Res(chromosome_index)] = ...
            calculate_fitness(X_train, Y_train, X_test, Y_test, regression_method, parameters, selected_features, dataset_name);
    end
    %%%%% sort chromosomes based on their fitness:
    [fitness_sorted_chromosomes, indices_sorted_chromosomes] = sort(fitness, 'ascend');
    if generation == 1 || (fitness_sorted_chromosomes(1) < best_fitness)
        best_fitness = fitness_sorted_chromosomes(1);
        best_chromosome_index = indices_sorted_chromosomes(1);
        best_chromosome = chromosomes(best_chromosome_index, :);
        best_selected_features = best_chromosome(end-number_of_features+1 : end);
        best_number_of_removed_features = sum(best_selected_features == 0);
        best_MMRE = MMRE(best_chromosome_index);
        best_PRED_25 = PRED_25(best_chromosome_index);
        best_Sum_Ab_Res = Sum_Ab_Res(best_chromosome_index);
        best_Med_Ab_Res = Med_Ab_Res(best_chromosome_index);
        best_SD_Ab_Res = SD_Ab_Res(best_chromosome_index);
    end
    %%%%% elitism replacement:
    population_of_elitism = elitism_population;
    indices_of_elitism = indices_sorted_chromosomes(1:population_of_elitism);
    chromosomes_next_generation = chromosomes(indices_of_elitism, :);
    chromosomes(indices_of_elitism, :) = [];
    fitness(indices_of_elitism) = [];
    probability_crossover = (100 ./ fitness) / sum(100 ./ fitness);
    cumulative_probability_crossover = cumsum(probability_crossover);
    for mating_index = 1:size(chromosomes, 1)/2
        %%%%% roulette wheel selection of parents:
        while true
            p = (rand - cumulative_probability_crossover);
            p = find(sign(p) == -1);
            index_of_parent1 = p(1);
            p = (rand - cumulative_probability_crossover);
            p = find(sign(p) == -1);
            index_of_parent2 = p(1);
            if index_of_parent1 ~= index_of_parent2
                break
            end
        end
        parent1 = chromosomes(index_of_parent1, :);
        parent2 = chromosomes(index_of_parent2, :);
        %%%%% two point random crossover:
        if rand <= crossover_rate
            cross_points = round(1 + ((dimension-1) - 1)*rand(2,1));
            cross_points = sort(cross_points, 'ascend');
            child1 = [parent1(1:cross_points(1)), parent2(cross_points(1)+1:cross_points(2)), parent1(cross_points(2)+1:end)];
            child2 = [parent2(1:cross_points(1)), parent1(cross_points(1)+1:cross_points(2)), parent2(cross_points(2)+1:end)];    
        else
            child1 = parent1;
            child2 = parent2;
        end
        %%%%% mutation:
        if rand <= mutation_rate
            child1 = ~child1;
        end
        if rand <= mutation_rate
            child2 = ~child2;
        end
        %%%%% add children to the next generation:
        chromosomes_next_generation(end+1, :) = child1;
        chromosomes_next_generation(end+1, :) = child2;
    end
    %%%%% update chromosomes:
    chromosomes = chromosomes_next_generation;
    %%%%% mutate invalid chromosomes (if all features are removed):
    for chromosome_index = 1:population_size  
        selected_features = chromosomes(chromosome_index, end-number_of_features+1 : end);
        if sum(selected_features) == 0
            a_feature_to_be_mutated = round(1 + (number_of_features - 1) * rand);
            chromosomes(chromosome_index, end-number_of_features+a_feature_to_be_mutated) = 1;
        end
    end
end

end