function [] = SimulatedAnnealing(Prob, N, prescription, SaveType)

%
% We incoroporate several variants of sa.
%
%	1 - as described in Djajaputra, Wu, Wu, and Mohan and in 
%		Stein, Mohan, Wang, Borfeld, Wu, Preiser, Ling and Schlegel
%		This tech. adjusts each angle simultaneously, which keeps
%		the spacing of the original set.
%	2 - adjust a single angle at a time, decrease the cooling structure
%		after N iterations.
%	3 - adjust all angles simultaneously but with different
%		variations.
%

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

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

%
% Cycle through the SA techniques
% 
for SAtype = 1:3

if SAtype == 1

   %grab system clock
   t0 = clock;

   %Initiate the starting angles as evenly spaced, rooted at 0
   Angles = round(1:Prob.NumAngles/N(n):Prob.NumAngles);

   % Set the Initial Temperature
   InitTemp = 1000;
   T = InitTemp;

   % Get an initial evaluation of the judgment function
   TmpProb = RestrictedSolver(Prob, prescription, Angles);

   % Stop when improvement has not been observed in the last N(n) updates.
   LastUpdate = 0;
   TempCount = 1;
   while LastUpdate < 100
      % Initate temporary angles
      TmpAngles = Angles;
   
      % Check the value of T
      if T < 10^(-10)
          T2 = 10^(-10);
      else
          T2 = T;
      end
      % randomly select an adjustment
      Adjustment = round(randraw('cauchy',[0, T2], 1));
   
      % Make the Update
      TmpAngles = mod(TmpAngles + Adjustment, Prob.NumAngles);
      TmpAngles(find(TmpAngles == 0)) = Prob.NumAngles;
   
      % Lets see if we accept the move.
      CandidateTmpProb = RestrictedSolver(Prob, prescription, TmpAngles);
      if CandidateTmpProb.fval < TmpProb.fval
         Angles = TmpAngles;
         TmpProb = CandidateTmpProb;
         LastUpdate = 0;
      elseif CandidateTmpProb.fval >= TmpProb.fval
         % grab the probability of acceptance.
         ProbAccept = randraw('exp', 1/T, 1);
         if ProbAccept < T*exp(-1 / T)
            Angles = TmpAngles;
            TmpProb = CandidateTmpProb;
            LastUpdate = 0;
         else
            LastUpdate = LastUpdate + 1;
         end
       else
          LastUpdate = LastUpdate + 1;
       end

       % Update counter and Temp
       TempCount = TempCount + 1;
       T = T / (1 + TempCount/10000);
   end
   % Make some figures
   Figure(TmpProb, prescription, sprintf('SA1_NumBeams%d', N(n)),'SA', SaveType)

   % Write to the table
   LaTex('AddRow', 'Simulated Annealing 1', 'N/A', Angles, etime(clock,t0)/60 )

elseif SAtype == 2
   % grab system clock
   t0 = clock;
  
   %Initiate the starting angles as evenly spaced, rooted at 0
   Angles = round(1:Prob.NumAngles/N(n):Prob.NumAngles);

   % Set the Initial Temperature
   InitTemp = 1000;
   T = InitTemp;

   % Get an initial evaluation of the judgment function
   TmpProb = RestrictedSolver(Prob, prescription, Angles);

   % Stop when improvement has not been observed in the last N(n) updates.
   LastUpdate = 0;
   TempCount = 1;
   while LastUpdate < 100
      % Initate temporary angles
      TmpAngles = Angles;
   
      % Select N(n) angles to adjust
      for a = 1:N(n)
         % Select an angle to adjust
         AdjustIndex = round((N(n)-1)*rand(1,1)) + 1;
      
         % Check the value of T
         if T < 10^(-10)
            T2 = 10^(-10);
         else
            T2 = T;
         end
         % randomly select an adjustment
         Adjustment = round(randraw('cauchy',[0, T2], 1));
   
         % Make the Update
         TmpAngles(AdjustIndex) = ...
            mod(TmpAngles(AdjustIndex) + Adjustment, Prob.NumAngles);
         TmpAngles(find(TmpAngles == 0)) = Prob.NumAngles;
      end
   
      % Lets see if we accept the move.
      CandidateTmpProb = RestrictedSolver(Prob, prescription, TmpAngles);
      if CandidateTmpProb.fval < TmpProb.fval
         Angles = TmpAngles;
         TmpProb = CandidateTmpProb;
         LastUpdate = 0;
      elseif CandidateTmpProb.fval >= TmpProb.fval
         % grab the probability of acceptance.
         ProbAccept = randraw('exp', 1/T, 1);
         if ProbAccept < T*exp(-1 / T)
            Angles = TmpAngles;
            TmpProb = CandidateTmpProb;
            LastUpdate = 0;
         else
            LastUpdate = LastUpdate + 1;
         end
       else
          LastUpdate = LastUpdate + 1;
       end

       % Update counter and Temp
       TempCount = TempCount + 1;
       T = T / (1 + TempCount/10000);
   end 
   % Make some figures
   Figure(TmpProb, prescription, sprintf('SA2_NumBeams%d', N(n)),'SA', SaveType)
   % Write to the table
   LaTex('AddRow', 'Simulated Annealing 2', 'N/A', Angles, etime(clock, t0)/60 )

elseif SAtype == 3
   % grab system clock
   t0 = clock;

   %Initiate the starting angles as evenly spaced, rooted at 0
   Angles = round(1:Prob.NumAngles/N(n):Prob.NumAngles);

   % Set the Initial Temperature
   InitTemp = 1000;
   T = InitTemp;

   % Get an initial evaluation of the judgment function
   TmpProb = RestrictedSolver(Prob, prescription, Angles);

   % Stop when improvement has not been observed in the last 100*N(n) updates.
   LastUpdate = 0;
   TempCount = 1;
   while LastUpdate < 100
      % Initate temporary angles
      TmpAngles = Angles;
    
      % Check the value of T
      if T < 10^(-10)
         T2 = 10^(-10);
      else
         T2 = T;
      end
      % randomly select an adjustment
      Adjustment = round(randraw('cauchy',[0, T2], length(TmpAngles) ));
   
      % Make the Update
      TmpAngles = mod(TmpAngles + Adjustment', Prob.NumAngles);
      TmpAngles(find(TmpAngles == 0)) = Prob.NumAngles;
   
      % Lets see if we accept the move.
      CandidateTmpProb = RestrictedSolver(Prob, prescription, TmpAngles);
      if CandidateTmpProb.fval < TmpProb.fval
         Angles = TmpAngles;
         TmpProb = CandidateTmpProb;
         LastUpdate = 0;
      elseif CandidateTmpProb.fval >= TmpProb.fval
         % grab the probability of acceptance.
         ProbAccept = randraw('exp', 1/T, 1);
         if ProbAccept < T*exp(-1 / T)
            Angles = TmpAngles;
            TmpProb = CandidateTmpProb;
            LastUpdate = 0;
         else
            LastUpdate = LastUpdate + 1;
         end
       else
          LastUpdate = LastUpdate + 1;
       end

       % Update counter and Temp
       TempCount = TempCount + 1;
       T = T / (1 + TempCount/10000);
   end 
   % Make some figures
   Figure(TmpProb, prescription, sprintf('SA3_NumBeams%d', N(n)),'SA', SaveType)
   % Write to the table
   LaTex('AddRow', 'Simulated Annealing 3', 'N/A', Angles, etime(clock, t0)/60 )

else
   fprintf('RAD -> SAtype is not known \n')
end

% End for SAtype
end
% End for N(n)
end
