%Programmed by:         Dario Izzo (Advanced Concepts Team) 
%Date:                  03/2008
%Revision:              1 
%Tested by:             ----------
%
%This file implements a Simulated Annealing Global Optimiser
%
%Usage:  [f,x,nf,s] = SA_Bound(fname,x0,lb,ub,nfmax,y)
%
%Inputs:
%           fname   = function name to be optimised. this needs to be in the
%                       Matlab path and to be in the form J = myfun(x,y)
%                       where x is the decision vector and y parameter 
%           x0      = Initial point. Needs to be a column vector of
%                       dimension D 
%           lb      = lower bounds. Needs to be a column vector of
%                       dimension D
%           ub      = upper bounds. Needs to be a column vector of
%                       dimension D
%           nfmax   = number of iterations
%           T0,Ti,Tf  = The annealing schedule goes from T0 to 0. 
%                       We simulate, of the whole annealing schedule,
%                       only from Ti to Tf
%           jump    =   This is the variation of the objective
%                       function that results, at T0 in a 36% acceptance rate
%           y       = parameter to be passed to fname
%
%Example:
%           [f,x,nf,s] = SA_Bound('mga',x0,lb,ub,20000,100,100,0,10,MGAproblem)
%
%this simulates the whole annealing procedure in 20000 iteration from 100
%to 0 degrees. Later we can reanneal from 10 to 0 using
%           [f,x,nf,s] = SA_Bound('mga',x0,lb,ub,20000,100,10,0,10,MGAproblem)

function [f,x,nf,s] = SA_Bound(fname,x0,lb,ub,nfmax,T0,Ti,Tf,jump,y,sb,eb)

%Some optimiser option
refresh = 200;              %write output every 200 function evaluation
D = length(x0);             %problem dimension
T = Ti;                     %current temperature
 %This is the variation of the objective
                            %   function that results, at high temperatures
                            %   in a 36% acceptance rate


% Initial state, energy.
s = x0;
e = feval(fname,x0,y);
nf = 0;                     % Energy evaluation count (will be increased by one at the end)


% Initial "best" solution is x0 if sb and eb are not passed to the function
% (this allows to call the SA_Bound in cooperative solvers)
if nargin < 11
sb = s;
eb = e;
end

%Main algorithm cycle
while nf < nfmax
    
    sn = neighbour(s,lb,ub);    % Pick some neighbour within the bounds lb ub
    en = feval(fname,sn,y);     % Compute its energy.
    if en < eb                  % Is this a new best?                    
        sb = sn;                % Yes, save it.
        eb = en;
    end  
    P = Prob(e, en, T);         %Metropolis criteria
    if P > rand                 %Should we move to it?
        %Yes, change state.
        s = sn;
        e = en;
    end
    % One more evaluation done
    nf =nf+1;
    T = (Ti-Tf)*(1-nf/nfmax) + Tf;    %Linear annealing schedule: will bring the
                                      %temperature from T0 to Tf in nfmax
   
	
    %print result every refresh iterations
    if (rem(nf,refresh) == 0)
       fprintf(1,'Iteration: %d,  Best: %f, Temperature: %f, Current %f\n',nf,eb,T,e);
       for n=1:D
           fprintf(1,'best(%d) = %f\n',n,sb(n));
       end
    end
end %Main Cycle
f=eb;
x=sb;
nf = nf+1;
    
    %This function evaluates the neighbour of the current point, it
    %implements a variable size neighbour so that at T0 the entire space is
    %sampled, at small T only a small neighbourhood is sampled
    function sn = neighbour(s,lb,ub)
        coeff = T/T0 * (ub-lb);
        ds = (0.5-rand(D,1))*2;
        ds = coeff.*ds'; 
        sn = s + ds;
        
        %bound correction
        for j = 1:D
             if (sn(j) > ub(j)) || (sn(j) < lb(j))
                    sn(j) = lb(j) + rand * (ub(j) - lb(j));
             end
        end
        
    end

    %metropolis function
    function P=Prob(e,en,T)
        P = exp(-(en-e)/T*T0/jump);
        if isnan(P)
            P=0
        end
    end
    
end