% Laboratory 8: Population Genetics and Evolution
% thanks to Tim Towler for giving me the problem statement
% coded by Kris Stewart, 6/18/93
%
% this version - Genetics_lessmem
%
% has a reduced memory requirement because we do not
% store the history of the population, only the final
% outcome.  Therefore only need 2 vectors of length
% popsize for the entire simulation.
%
clear
disp('genetic1.m = updated version - uses rand <= 0.5 both tests ');
disp('and indexed population count - added timer');
popsize=input('how many individuals in population (even number)?');
numgen=input('how many generations to you want to run through?');
shuffle=input('enter 0 for shuffle each generation; else 1');
%
% no shuffling is not a good idea
%
% there will be 3 gene types. p(i)=0 for aa
%                                 =1 for aA
%                                 =2 for AA
% with the distribution 1:2:1 (i.e. aA and Aa are the same type)
%
% shuffle the indices into the population data
for i=1:popsize,
    index(i)=i;
end

% rand('dist'); % current distribution in a string
% initially all of type aA
p = ones(popsize,1);
% history = p;
for k=1:numgen,

titl = ['Genotype Data (genetic1); popsize=',num2str(popsize),'; Generation: ',num2str(k)];
hist(p),title(titl),pause(1)


    for i = 1:2:popsize,
% assume mating occurs with "the girl next door" after shuffling
       for j = 1:2,
% each pair produces 2 offspring
%
% each individual donates one of their genes, chosen at random
% the a gene is denoted as 0; the A gene as 1
           if p(index(i)) == 0,
              gene1(j) = 0;
           elseif p(index(i)) == 2,
              gene1(j) = 1;
           else
% randomly select which gene to contribute to offspring
% rand with no arguments is a scalar uniform in (0,1)
% the logical value return by the test rand > 0.5 corresponds
% to the gene selected 0 = false = a; 1 = true = A
%
% will try to balance out the inequality test by biasing the
% second gene selection the other way
              gene1(j) = rand <= 0.5;
           end
% repeat for "the girl next door"
           if p(index(i+1)) == 0,
              gene2(j) = 0;
           elseif p(index(i+1)) == 2,
              gene2(j) = 1;
           else
% randomly select which gene to contribute to offspring
              gene2(j) = rand <= 0.5;
           end
      end % loop on j for the 2 individuals in each pair
%
% we now know the gene to be contributed by each parent to
% produce 2 offspring.  the model assumes a very short
% reproductive cycle, i.e. the parents will not reproduce
% again, there the 2 offspring will replace the 2 parents
% in the gene pool
%
% gene for first child
    p(index(i)) = gene1(1) + gene2(1);
% gene for second child
    p(index(i+1)) = gene1(2) + gene2(2);
  end  % loop on i for each individual in population
%
% plot the distribution of genes at each generation
%
%plot(p), pause
% history = [history p];

if shuffle == 0
   i=1;
   for x=popsize:-1:1,
       num = ceil(rand*x);
       newindex(i) = index(num);
       i = i+1;
       index(num) = index(x);
   end % for x=popsize:-1:1
% update index array to new, randomized values
index = newindex;
end % if shuffle block

end % loop on k for each generation (numgen)
%plot(history)
% got to be a slicker way to do this - but sum up the
% diversity of the final generation
count0=0; count1=0; count2=0;
for i=1:popsize,
    if p(i)==0,
       count0 = count0+1;
    elseif p(i)==1,
       count1 = count1+1;
    else
       count2 = count2+1;
    end
end
if shuffle == 0,
    disp('With shuffling')
else
    disp('No shuffling')
end
% this was my first attempt at getting i/o with labels so that
% the results could be easily read.  it worked and i encourage
% you, especially on the early tries, just get the values out
%
%disp(' For a population of '), disp(popsize), disp(' individuals')
%disp(' after '), disp(numgen), disp(' generations')
%disp(' We end up with')
%disp(count0), disp( ' of type aa ' )
%disp(count1), disp( ' of type Aa   and ')
%disp(count2), disp( ' of type AA - was this pop large enough?')
%
% but now i want them to be pretty - so mucking around with string
%
str1=['For a population of ' num2str(popsize) ' individuals'];
str2=['   after ' num2str(numgen) ' generations '];
str3=['We have :'];
str4=['           ' num2str(count0) ' of type aa'];
str5=['           ' num2str(count1) ' of type aA'];
str6=['           ' num2str(count2) ' of type AA'];
str7=['Do you think the population was large enough?'];
str8=['Do you think the simulation ran for enough generations?'];

disp(str1), disp(str2), disp(str3), disp(str4)
disp(str5), disp(str6), disp(str7), disp(str8)

if shuffle == 0
titl = ['Genotype Data (genetic1), with shuffling; popsize=',num2str(popsize),'; numgen=',num2str(numgen)];
else
titl = ['Genotype Data (genetic1), no shuffling; popsize=',num2str(popsize),'; numgen=',num2str(numgen)];
end

hist(p),title(titl)