function [xmin,fmin,nmin,xnow,fnow] = sim_anneal(myf,a,b,T,alpha,epsT)

% function [xmin,fmin,nmin] = sim_anneal(myf,a,b,T,alpha,epsT)
%
% This function minimizes the function myf (of a single variable)
% over the interval [a,b]
% using simulated annealing.
%
%
% Inputs:
%   myf    a "handle" for a Matlab function
%   a      left  endpoint of the interval
%   b      right endpoint of the interval
%   T      initial "temperature" for annealing
%   alpha  reduction factor for temperature at each iteration
%   epsT   final value of T
%
%  Outputs:
%    xmin  the computed minimizer
%    fmin  the computed minimum myf(xmin)
%    xnow  the final point tested in the algorithm
%    fnow  the function value at this point myf(xnow)
%    nmin  the number of function evaluations used
%
%  sim_anneal Dianne P. O'Leary 10/2006


% Initialize our best guess at the minimizer to be the left endpoint.

xnow = a;
fnow = feval(myf,a);
xmin = a;
fmin = fnow;

% Initialize an array of Nsamp random starting points.

x0 =  rand(1000,1);
jsamp = 0;
nmin  = 0;
factor = 1;

% Until the temperature is at "freezing",

while (T > epsT)

% Generate a random step from the current x (xnow).
% Make sure that the new point satisfies the bounds.
% The size of the random step decreases with T.

   jsamp = jsamp + 1;
   if (jsamp > 998)
      x0 = rand(1000,1);
      jsamp = 1;
   end
   x = xnow + factor*(x0(jsamp)-.5);
   x = min(b,x);
   x = max(a,x);
   nmin  = nmin + 1;

   fx = feval(myf,x);
   t = fx - fnow;

% If the new point has a lower function value, replace xnow by it.
% Otherwise, replace it with a probability depending on how
% much the function increases, and the current temperature.

   if (t < 0)
      xnow = x;
      fnow = fx;
   else
      jsamp = jsamp + 1;
      if (x0(jsamp) < exp(-t/T))
         xnow = x;
         fnow = fx;
      end
   end
   T = alpha * T;
   factor = sqrt(alpha) * factor;

% See if x is the best point seen yet.

   if (fx < fmin)
      xmin = x;
      fmin = fx;
   end

% For debugging:
%     disp([x,fx,xnow,fnow,xmin,fmin,factor])
end

