anneal(prob)

function r = anneal(prob)
% function [minimum,fval] = anneal(loss, parent, options)
% ANNEAL  Minimizes a function with the method of simulated annealing
% (Kirkpatrick et al., 1983)
%
%  ANNEAL takes three input parameters, in this order:
%
%  LOSS is a function handle (anonymous function or inline) with a loss
%  function, which may be of any type, and needn't be continuous. It does,
%  however, need to return a single value.
%
%  PARENT is a vector with initial guess parameters. You must input an
%  initial guess.
%
%  OPTIONS is a structure with settings for the simulated annealing. If no
%  OPTIONS structure is provided, ANNEAL uses a default structure. OPTIONS
%  can contain any or all of the following fields (missing fields are
%  filled with default values):
%
%       Verbosity: Controls output to the screen.
%                  0 suppresses all output
%                  1 gives final report only [default]
%                  2 gives temperature changes and final report 
%       Generator: Generates a new solution from an old one.
%                  Any function handle that takes a solution as input and
%                  gives a valid solution (i.e. some point in the solution
%                  space) as output.
%                  The default function generates a row vector which
%                  slightly differs from the input vector in one element:
%                  @(x) (x+(randperm(length(x))==length(x))*randn/100)
%                  Other examples of possible solution generators:
%                  @(x) (rand(3,1)): Picks a random point in the unit cube
%                  @(x) (ceil([9 5].*rand(2,1))): Picks a point in a 9-by-5
%                                                 discrete grid
%        InitTemp: The initial temperature, can be any positive number.
%                  Default is 1.
%        StopTemp: Temperature at which to stop, can be any positive number
%                  smaller than InitTemp. 
%                  Default is 1e-8.
%         StopVal: Value at which to stop immediately, can be any output of
%                  LOSS that is sufficiently low for you.
%                  Default is -Inf.
%       CoolSched: Generates a new temperature from the previous one.
%                  Any function handle that takes a scalar as input and
%                  returns a smaller but positive scalar as output. 
%                  Default is @(T) (.8*T)
%      MaxConsRej: Maximum number of consecutive rejections, can be any
%                  positive number.
%                  Default is 1000.
%        MaxTries: Maximum number of tries within one temperature, can be
%                  any positive number.
%                  Default is 300.
%      MaxSuccess: Maximum number of successes within one temperature, can
%                  be any positive number.
%                  Default is 20.
%
%
%  Usage:
%     [MINIMUM,FVAL] = ANNEAL(LOSS,NEWSOL,[OPTIONS]);
%          MINIMUM is the solution which generated the smallest encountered
%          value when input into LOSS.
%          FVAL is the value of the LOSS function evaluated at MINIMUM.
%     OPTIONS = ANNEAL();
%          OPTIONS is the default options structure.
%
%
%  Example:
%     The so-called "six-hump camelback" function has several local minima
%     in the range -3<=x<=3 and -2<=y<=2. It has two global minima, namely
%     f(-0.0898,0.7126) = f(0.0898,-0.7126) = -1.0316. We can define and
%     minimise it as follows:
%          camel = @(x,y)(4-2.1*x.^2+x.^4/3).*x.^2+x.*y+4*(y.^2-1).*y.^2;
%          loss = @(p)camel(p(1),p(2));
%          [x f] = ANNEAL(loss,[0 0])
%     We get output:
%               Initial temperature:     	1
%               Final temperature:       	3.21388e-007
%               Consecutive rejections:  	1027
%               Number of function calls:	6220
%               Total final loss:        	-1.03163
%               x =
%                  -0.0899    0.7127
%               f =
%                  -1.0316
%     Which reasonably approximates the analytical global minimum (note
%     that due to randomness, your results will likely not be exactly the
%     same).

%  Reference:
%    Kirkpatrick, S., Gelatt, C.D., & Vecchi, M.P. (1983). Optimization by
%    Simulated Annealing. _Science, 220_, 671-680.
%   Joachim Vandekerckhove
%   [e-mail]
%   $Revision: v5 $  $Date: 2006/04/26 12:54:04 $
r = oor('authors: Joachim Vandekerckhove [e-mail]', ...
    'alg:Kirkpatrick, S., Gelatt, C.D., & Vecchi, M.P. (1983). Optimization by Simulated Annealing',...
    'lastChanger: Dmitrey');

% neighborhood space function
% generator = @(x)(x+transpose(randperm(length(x))==length(x))*randn/100);
% by Dmitrey: calling newsol as separate func is too slow
% @(x) (x+sparse(ceil(rand*length(x)), 1, randn/100, length(x), 1, 1);
% is the same but a little bit faster
% I implenent it directly

Tinit = 1;        % initial temp
minT =1e-8;         % stopping temp
cool = @(T) (.8*T);        % annealing schedule
max_consec_rejections = 1000;
max_try = 500;
max_success = 20;
k = 1;                           % boltzmann constant
%checking are any solver parameters is structure prob.(name of solver)
ExtractRoutineParamsFromProb;%script-file
%if present - they will replace defaults 

userGeneratorFlag = 0;
if exist('generator','var'); userGeneratorFlag=1; end


% counters etc
itry = 0;
success = 0;
consec = 0;
T = Tinit;
r.xf = prob.x0;
r.x = prob.x0;
initenergy = prob.f(r.xf, prob);
r.ff = initenergy;
r.f = initenergy;
total = 0;
while 1
    itry = itry+1; % just an iteration counter
    current = r.x; 
    
    % % Stop / decrement T criteria
    if itry >= max_try || success >= max_success
        r = prob.iterfcn(prob, r);
        if r.istop; return; end
        if consec >= max_consec_rejections
            r.istop = 12; return 
        elseif T < minT 
            r.istop = 9;return
        else
            T = cool(T);  % decrease T according to cooling schedule
            total = total + itry;
            itry = 1;
            success = 1;
        end
    end
    
    %newparam = newsol(current);
    % by Dmitrey:
    if userGeneratorFlag
        newparam = generator(current);
    else
        %this inline is much more faster than original anneal.m version
        newparam = current;
        rand_ind = ceil(rand*length(current));
        newparam(rand_ind) = newparam(rand_ind)+randn/100;
    end

    
    newenergy = prob.f(newparam, prob);
   
    if r.f-newenergy > min(prob.TolFun/2, 1e-6)
        r.x = newparam;
        r.f = newenergy;
        success = success+1;
        consec = 0;
    else
        if rand < exp( (r.f-newenergy)/(k*T) )
            r.x = newparam;
            r.f = newenergy;
            success = success+1;
        else
            consec = consec+1;
        end
    end
    if r.f<r.ff
        r.ff = r.f;
        r.xf = r.x;
    end
end