function [Soln] = GA(Prob, N, prescription, SaveType, Soln)


%
%  Two different genetic algorithsm.
%
%	1 - as described by Haas, Burnham, and Mills - An individual
%		is a collection of anlges.
%		collection of angles
%	2 - Hpu, Wang, and Chen - An individual is a collection of angles
%

% Set Update Limit
UpdateLimit = 10;

%
% Make the Dir. for images
%
if ~isdir('../images/GA')
   !mkdir ../images/GA
end

%
% Cycle through the different angle numbers
%
for n = 1:length(N)

%
% Cycle through the SA techniques
%
for GAtype = 1:2

if GAtype == 1
   % grab system clock
   t0 = clock;
   % Randomly generate angle sets
   SizeOfPopulation = round(Prob.NumAngles/(2*N(n)));
   Population = zeros(SizeOfPopulation,N(n)+1);
   for p = 1:SizeOfPopulation
      % Generate an individual
      Population(p,2:N(n)+1) = ...
            mod(round((Prob.NumAngles-1)*rand(1,N(n)) + 1), Prob.NumAngles);
      Population(p,find(Population(p,:) == 0)) = Prob.NumAngles;
      % Grab the judgement function
      TmpProb = RestrictedSolver(Prob, prescription, Population(p,2:N(n)+1));
      % Update the inidividuals judgement value
      Population(p,1) = TmpProb.fval;
   end

   % Rank the individuals
   [FValues, FIndex] = sort(Population(:,1));
   Population = Population(FIndex, :);
   BestValue = Population(1,1);
   BestAngles = Population(1,2:N(n)+1);

   % set the number of operations to undergo
   NumberOfOperations = round(SizeOfPopulation / 4);

   %Stop when no improvement has been found in the last 100 updates
   LastUpdate = 0;
   IterationCount = 0;
   while LastUpdate < UpdateLimit
   IterationCount = IterationCount + 1;
   % Assume that no update is made
   Improved = 0;

   %Randomly select operations and populations
   OperationType = round(4*rand(NumberOfOperations,1))+1;
   PopsToAdjust = round((SizeOfPopulation-1)*rand(NumberOfOperations,1))+1;

   for d = 1:NumberOfOperations

   %Perform the operation
   switch OperationType(d)
      case 1
         % Generate a directions set
         IndPrime = round((SizeOfPopulation-1)*rand(1,1)+1);
         DirectionAngles = Population(IndPrime, 2:N(n)+1);
         % Get a step length
         Step = 1.5*rand(1,1) - 0.25;
         % Update Angles
         Angles = Population(PopsToAdjust(d), 2:N(n)+1);
         Angles = round(mod(Angles + Step * (DirectionAngles - Angles), ...
	    Prob.NumAngles));
         Angles(find(Angles == 0)) = Prob.NumAngles;
         % Grab the judgement function
         CandidateTmpProb = RestrictedSolver(Prob, prescription, Angles);
         % Put Everything back into the population
         Population(PopsToAdjust(d),:) = [CandidateTmpProb.fval, Angles];
         % Check if angle collection is better
         if CandidateTmpProb.fval < BestValue
            BestAngles = Angles;
            BestValue = CandidateTmpProb.fval;
            TmpProb = CandidateTmpProb;
            Improved = 1;
         end
      case 2
         Angles = Population(PopsToAdjust(d), 2:N(n)+1);
         % Select an Angle and shift
         SelectedIndex = round((N(n)-1)*rand(1,1)) + 1;
         Shift = round(6*rand(1,1) - 3);
         % Update the Angle set
         Angles(SelectedIndex) = round(mod(Angles(SelectedIndex) + Shift, ...
                                     Prob.NumAngles));
         Angles(find(Angles == 0)) = Prob.NumAngles;
         % Grab the judgement function
         CandidateTmpProb = RestrictedSolver(Prob, prescription, Angles);
         % Put Everything back into the population
         Population(PopsToAdjust(d),:) = [CandidateTmpProb.fval, Angles];
         % Check if angle collection is better
         if CandidateTmpProb.fval < BestValue
            BestAngles = Angles;
            BestValue = CandidateTmpProb.fval;
            TmpProb = CandidateTmpProb;
            Improved = 1;
         end
      case 3
         Angles = Population(PopsToAdjust(d), 2:N(n)+1);
         % Select an Angle and shift
         SelectedIndex = round((N(n)-1)*rand(1,1)) + 1;
         Shift = round((Prob.NumAngles - 1)*rand(1,1)) + 1;
         % Update the Angle set
         Angles(SelectedIndex) = round(mod(Angles(SelectedIndex) + Shift, ...
                                     Prob.NumAngles));
         Angles(find(Angles == 0)) = Prob.NumAngles;
         % Grab the judgement function
         CandidateTmpProb = RestrictedSolver(Prob, prescription, Angles);
         % Put Everything back into the population
         Population(PopsToAdjust(d),:) = [CandidateTmpProb.fval, Angles];
         % Check if angle collection is better
         if CandidateTmpProb.fval < BestValue
            BestAngles = Angles;
            BestValue = CandidateTmpProb.fval;
            TmpProb = CandidateTmpProb;
            Improved = 1;
         end
      case 4
         Angles = Population(PopsToAdjust(d), 2:N(n)+1);
         % Select an Angle and shift
         SelectedIndex = round((N(n)-1)*rand(1,1)) + 1;
         % Update the Angle set
         Angles(SelectedIndex) = round(mod(Angles(SelectedIndex) + ...
                                 round(Prob.NumAngles/2), Prob.NumAngles));
         Angles(find(Angles == 0)) = Prob.NumAngles;
         % Grab the judgement function
         CandidateTmpProb = RestrictedSolver(Prob, prescription, Angles);
         % Put Everything back into the population
         Population(PopsToAdjust(d),:) = [CandidateTmpProb.fval, Angles];
         % Check if angle collection is better
         if CandidateTmpProb.fval < BestValue
            BestAngles = Angles;
            BestValue = CandidateTmpProb.fval;
            TmpProb = CandidateTmpProb;
            Improved = 1;
         end
      case 5
         Angles = Population(PopsToAdjust(d), 2:N(n)+1);
         % Select an Angle and shift
         SelectedIndex = round((N(n)-1)*rand(1,1)) + 1;
         % Select a uniform distribution
         DirectionAngles = round(1:Prob.NumAngles/N(n):Prob.NumAngles);
         % Update Angles
         Angles = round(mod(Angles(SelectedIndex) +  DirectionAngles, ...
                      Prob.NumAngles));
         Angles(find(Angles == 0)) = Prob.NumAngles;
         % Grab the judgement function
         CandidateTmpProb = RestrictedSolver(Prob, prescription, Angles);
         % Put Everything back into the population
         Population(PopsToAdjust(d),:) = [CandidateTmpProb.fval, Angles];
         % Check if angle collection is better
         if CandidateTmpProb.fval < BestValue
            BestAngles = Angles;
            BestValue = CandidateTmpProb.fval;
            TmpProb = CandidateTmpProb;
            Improved = 1;
         end
      otherwise
         LogMesg('RAD -> No Genetic type 1 move available')
      end

      % end for loop for number of operations
      end
      
      % Check to see if an update was made
      if Improved
         LastUpdate = 0;
      else
         LastUpdate = LastUpdate + 1;
      end
      % End While LastUpdate loop
      end

% Catalog some times
eval(sprintf('Soln.BeamSelector.GA1.B%d.SelectionTime = %6.4f;', ...
   N(n), etime(clock,t0) ));
t0 = clock;
TmpProb = RestrictedSolver(Prob, prescription, BestAngles);
eval(sprintf('Soln.BeamSelector.GA1.B%d.SolutionTime = %6.4f;', ...
   N(n), etime(clock,t0) ));

% Balance the Solution
t0 = clock;
TmpProb.ThePlan = BalancedPDF(TmpProb);
eval(sprintf('Soln.BeamSelector.GA1.B%d.BalancingTime = %6.4f;', ...
   N(n), etime(clock,t0) ));

% Make some figures
Figure(TmpProb, prescription, sprintf('GA1_NumBeams%d', N(n)),'GA', SaveType)

% Write to the table
LaTex('AddRow', 'Genetic Algorithm 1', 'N/A', BestAngles, etime(clock,t0) )

% Catalog some Soln Stuff
eval(sprintf('Soln.BeamSelector.GA1.B%d.Angles = BestAngles;', N(n) ));
eval(sprintf('Soln.BeamSelector.GA1.B%d.Iterations = IterationCount;', ...
   N(n) ));

elseif GAtype == 2
   % Grab system clock
   t0 = clock;
   % Randomly generate random angle sets
   SizeOfPopulation = round(Prob.NumAngles/(2*N(n)));
   Population = zeros(SizeOfPopulation,N(n)+1);
   for p = 1:SizeOfPopulation
      % Generate an individual
      Population(p,2:N(n)+1) = ...
            mod(round((Prob.NumAngles-1)*rand(1,N(n)) + 1), Prob.NumAngles);
      Population(p,find(Population(p,:) == 0)) = Prob.NumAngles;
      % Grab the judgement function
      TmpProb = RestrictedSolver(Prob, prescription, Population(p,2:N(n)+1));
      % Update the inidividuals judgement value
      Population(p,1) = TmpProb.fval;
   end

   % Rank the individuals
   [FValues, FIndex] = sort(Population(:,1));
   Population = Population(FIndex, :);
   BestValue = Population(1,1);
   BestAngles = Population(1,2:N(n)+1);

   % set the number of operations to undergo
   NumberOfOperations = round(SizeOfPopulation / 4);

   LastUpdate = 0;
   while LastUpdate < UpdateLimit

   % Make room for new individuals
   NewIndividuals = zeros(2*NumberOfOperations, N(n)+1);
   for i = 1:NumberOfOperations
      % Find out which two to play with
      SelectedIndividuals = round(randraw('tri', ...
         [1,  SizeOfPopulation - .001, SizeOfPopulation], 2));
      NewInd1 = Population(SelectedIndividuals(1),:);
      NewInd2 = Population(SelectedIndividuals(2),:);
      % Perform a crossover with probability 1/2
      if rand(1,1) < 0.5
         % Select where to cut
         CutIndex = round((N(n) - 1)*rand(1,1) + 2);
         NewInd1 = [Population(SelectedIndividuals(1),1:CutIndex), ...
                    Population(SelectedIndividuals(2),CutIndex+1:N(n)+1)];
         NewInd2 = [Population(SelectedIndividuals(2),1:CutIndex), ...
                    Population(SelectedIndividuals(1),CutIndex+1:N(n)+1)];
      end
      % Perform mutation with probability 1/2
      if rand(1,1) < 0.5
         % Randomly select a shift
         Shift = randraw('norm', [0, Prob.NumAngles/2], 1);
         % Select an Angle to shift
         SelectedIndividual = round((N(n) - 1)*rand(1,1) + 1);
         % Make the mutation
         NewInd1(SelectedIndividual) = ...
            round(mod(NewInd1(SelectedIndividual) + Shift, Prob.NumAngles));
         NewInd1(find(NewInd1 == 0)) = Prob.NumAngles;
      end
      if rand(1,1) < 0.5
         % Randomly select a shift
         Shift = randraw('norm', [0, Prob.NumAngles/2], 1);
         % Select an Angle to shift
         SelectedIndividual = round((N(n) - 1)*rand(1,1) + 2);
         % Make the mutation
         NewInd2(SelectedIndividual) = ...
            round(mod(NewInd2(SelectedIndividual) + Shift, Prob.NumAngles));
         NewInd2(find(NewInd2 == 0)) = Prob.NumAngles;
      end
      % Save the New Individuals
      TmpProb = RestrictedSolver(Prob, prescription, NewInd1(2:N(n)+1));
      NewIndividuals(2*i-1,:) = NewInd1;
      NewInd1(2*i-1,1) = TmpProb.fval;
      TmpProb = RestrictedSolver(Prob, prescription, NewInd2(2:N(n)+1));
      NewIndividuals(2*i,:) = NewInd2;
      NewInd1(2*i,1) = TmpProb.fval;
   end
   % Update the population and ReRank
   Population(SizeOfPopulation-(2*NumberOfOperations)+1: ...
              SizeOfPopulation,:) = NewIndividuals;
   [FValues, FIndex] = sort(Population(:,1));
   Population = Population(FIndex, :);

   %Check if we have improved
   if Population(1,1) < BestValue
      LastUpdate = 0;
      BestAngles = Population(1,2:N(n)+1);
      BestValue = Population(1,1);
   else
      LastUpdate = LastUpdate + 1;
   end
   % End While
   end

% Catalog some times
eval(sprintf('Soln.BeamSelector.GA2.B%d.SelectionTime = %6.4f;', ...
   N(n), etime(clock,t0) ));
t0 = clock;
TmpProb = RestrictedSolver(Prob, prescription, Angles);
eval(sprintf('Soln.BeamSelector.GA2.B%d.SolutionTime = %6.4f;', ...
   N(n), etime(clock,t0) ));

% Balance the Solution
t0 = clock;
TmpProb.ThePlan = BalancedPDF(TmpProb);
eval(sprintf('Soln.BeamSelector.GA1.B%d.BalancingTime = %6.4f;', ...
   N(n), etime(clock,t0) ));

% Make some figures
TmpProb = RestrictedSolver(Prob, prescription, BestAngles);
Figure(TmpProb, prescription, sprintf('GA2_NumBeams%d', N(n)),'GA', SaveType)

% Write to the table
LaTex('AddRow', 'Genetic Algorithm 2', 'N/A', BestAngles, etime(clock, t0) )

% Catalog some Soln Stuff
eval(sprintf('Soln.BeamSelector.GA2.B%d.Angles = BestAngles;', N(n) ));
eval(sprintf('Soln.BeamSelector.GA2.B%d.Iterations = IterationCount;', ...
   N(n) ));

else
   LogMesg('RAD -> Confused GAtype')
end

% End GAtype loop
end

% End N(n) loop
end
