function [P,best] = gen_alg(pop_size,chrom_len,chrom_vals,pm,pc,max_gen,disp_i)
    % /---------------------------------------------------------------------\
    % | Copyright (C) 2009  George Bezerra                                   |
    % |                                                                      |
    % | This program is free software: you can redistribute it and/or modify |
    % | it under the terms of the GNU General Public License as published by |
    % | the Free Software Foundation, either version 3 of the License, or    |
    % | (at your option) any later version.                                  |
    % |                                                                      |
    % | This program is distributed in the hope that it will be useful,      |
    % | but WITHOUT ANY WARRANTY; without even the implied warranty of       |
    % | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        |
    % | GNU General Public License for more details.                         |
    % |                                                                      |
    % | You should have received a copy of the GNU General Public License    |
    % | along with this program.  If not, see <http://www.gnu.org/licenses/>.|
    % \ ---------------------------------------------------------------------/ 
    %
    % Inputs:
    % pop_size   => population size
    % chrom_len  => chromosome length
    % chrom_vals => chromosone values (0,1)
    % pm         => probability of mutation
    % pc         => probability of crossover
    % max_gen    => maximum number of generations
    % disp_i     => Number of generations between output
    %
    %
    % Outputs:
    % P    => population
    % best => best individual of the population
    %
    % suggested run: [P,best] = ga(100,100,[0,1],0.01,0.5,100,10);


    % INITIALIZE POPULATION
    P = initialize(pop_size,chrom_len,chrom_vals);

    % EVALUATION
    fit = maxrocks_fitness(P);

    gen = 1;
    max_fit = zeros(max_gen,1);
    mean_fit = zeros(max_gen,1);
    high_95fit = zeros(max_gen,1);
    low_95fit = zeros(max_gen,1);
    while(gen<=max_gen)

        % SELECTION
        P = tournament_selection(P,fit,2);

        % CROSSOVER
        P = two_point_crossover(P,pc);

        % MUTATION
        P = point_mutation(P,pm,chrom_vals);

        % EVALUATION
        fit = maxrocks_fitness(P);

        % record data
        sfit = sort(fit,'descend');
        max_fit(gen) = sfit(1);
        mean_fit(gen) = mean(fit);
        high95_fit(gen) = sfit(round(.05*pop_size));
        low95_fit(gen) = sfit(round(.95*pop_size));

        % display information
        disp_info(gen,max_fit(gen),disp_i);

        gen = gen+1;
    end

    disp(sprintf('Generation: %d',gen));
    disp(sprintf('Best fitness: %d\n',max_fit(end)));

    % plot evolution curve
    plot_ci(1:length(mean_fit),low95_fit,high95_fit,[0,0,0],.1)
    hold on;
    plot(1:length(max_fit), max_fit,'b');
    plot(1:length(mean_fit), mean_fit,'g');

    hold off;
    xlabel('Generations');
    ylabel('Fitness');
    BL = legend('90% Fitness Confidence Interval','Best fitness','Average fitness','Location','SouthEast');
    PatchInLegend = findobj(BL, 'type', 'patch');
    set(PatchInLegend, 'facea', .1);

    % output best individual
    [m,ind] = max(fit);
    best = P(ind,:);
end

function [] = plot_ci(x,lower,upper,color,alpha)
    set(fill([x,x(end:-1:1)],[upper,lower(end:-1:1)],color),'EdgeColor',color,'FaceAlpha',alpha,'EdgeAlpha',alpha);
end

function [P] = initialize(pop_size,chrom_length,chrom_vals)
    P = chrom_vals(randi([1,length(chrom_vals)],[pop_size,chrom_length]));
end

%Rewrite this function, this only finds which one has the most Rocks!
function [fit] = maxrocks_fitness(P)
    for i=1:size(P,1)
        fit(i) = length(find(P(i,:)=='R'));
    end
end

function [P_new] = tournament_selection(P,fit,tourn_size)
    for i=1:size(P,1)
        % t gets indexes into the individuals chosen for a tournament
        t = randi([1,size(P,1)],[1,tourn_size]);
        [max_fit,winner] = max(fit(t));
        P_new(i,:) = P(t(winner),:);
    end
end

function [P_new] = two_point_crossover(P,pc)
    mating_list = randperm(size(P,1));
    P_new = [];
    while ~isempty(mating_list)
        % pair gets the first two individuals in mating list and they
        % are removed from mating_list
        pair = mating_list(1:2);
        mating_list(1:2) = [];
        if rand<pc
            % choose 2 crossover points
            crossover_points = randi([1,size(P,2)],[1,2]);
            point1 = min(crossover_points);
            point2 = max(crossover_points);
            individual1 = P(pair(1),:);
            individual2 = P(pair(2),:);
            % Do the swap
            individual1(point1:point2) = P(pair(2),point1:point2);
            individual2(point1:point2) = P(pair(1),point1:point2);
            P_new = [P_new;individual1;individual2];
        else
            P_new = [P_new;P(pair,:)];
        end
    end
end

function [P_new] = point_mutation(P,pm,chrom_vals)
    % true or false for whether we mutate stuff
    r = logical(random('Binomial',1,pm,size(P)));
    P_new = P;
    % note there is a 1/length(chrom_vals) we mutate to the same value
    P_new(r) = chrom_vals(randi([1,length(chrom_vals)],size(P(r),2)));
end

function [] = disp_info(gen,fit,step_size)
    if mod(gen,step_size)==0 
        disp(sprintf('Generation: %d',gen));
        disp(sprintf('Best fitness: %d\n',fit));
    end
end

function [P_new] = roulette_selection(P,fit)
    fit = (fit - min(fit)).^2;
    fit = cumsum(fit);
    fit = fit/max(fit);
    P_new = [];
    for i=1:size(P,1)
        f = find(fit>rand);
        P_new = [P_new;P(f(1),:)];
    end
end
