%       @mystate, mydomain, ...
%       @mycost, @myneighbor, ...
%       @metropolis, [], ...
%       @thermospeed, 0, ...
%       @hoffman, 0.75, 20, ...
%       @TinitWhite, [1.7, 3], ...
%       1, 0, 50, ...
%       0.3, 10, 0) ;
%
error(nargchk(21,21,nargin)) ;
%
% Check for valid input
%
if ~isa(newstate, 'function_handle')
    error('No function handle supplied for newstate') ;
end
classX = class(X) ;
sizeX = size(X) ;
if ~isa(cost, 'function_handle')
    error('No function handle supplied for cost') ;
end
if walkers < 1
    error('Number of walkers must be positive') ;
end
if ~isa(acceptrule, 'function_handle')
    error('No function handle supplied for acceptrule') ;
end
classQ = class(q) ;
sizeQ = size(q) ;
if ~isa(schedule, 'function_handle')
    error('No function handle supplied for schedule') ;
end
classP = class(P) ;
sizeP = size(P) ;
if isa(equilibrate, 'function_handle')
    hasEquilibrate = 1 ;
else
    hasEquilibrate = 0 ;
end
classC = class(C) ;
sizeC = size(C) ;
if maxsteps < 1
    error('maxsteps must be positive') ;
end
if isa(Tinit, 'function_handle')
    hasTinitMethod = 1 ;
else
    if ~isa(Tinit, 'numeric')
      error('No numeric value or function handle supplied for Tinit') ;
    end
    hasTinitMethod = 0 ;
end
if ~isa(r, 'numeric')
    error('No numeric value supplied for r') ;
end
if isa(Tfinal, 'function_handle')
    hasTfinalMethod = 1 ;
else
    if ~isa(Tfinal, 'numeric')
      error('No numeric value or function handle supplied for Tfinal') ;
    end
    hasTfinalMethod = 0 ;
end
if ~isa(f, 'numeric')
    error('No numeric value supplied for f') ;
end
if maxtemps < 1
    error('maxtemps must be positive') ;
end
if ~isa(v, 'numeric')
    error('No numeric value supplied for v') ;
end
if ~isa(bins, 'numeric')
    error('No numeric value supplied for bins') ;
end
if ~isa(e, 'numeric')
    error('No numeric value supplied for e') ;
end
%
%
if verbose
    newline = sprintf('\n') ;
    tab = sprintf('\t') ;
    [vnum, vdate] = satoolsversion ;
    disp(['SA Tools anneal. Version ', vnum, ', Last update ', vdate, '.', newline]) ;
    disp([tab, 'newstate = ', func2str(newstate)]) ;
    disp([tab, 'X is ', classX, ' of size ', num2str(sizeX)]) ;
    disp([tab, 'cost = ', func2str(cost)]) ;
    disp([tab, 'moveclass = ', func2str(moveclass)]) ;
    disp([tab, 'walkers = ', num2str(walkers)]) ;
    disp([tab, 'acceptrule = ', func2str(acceptrule)]) ;
    disp([tab, 'q = ', num2str(q)]) ;
    disp([tab, 'schedule = ', func2str(schedule)]) ;
    disp([tab, 'P = ', num2str(P)]) ;
    if hasEquilibrate
        disp([tab, 'equilibrate = ', func2str(equilibrate)]) ;
    else
        disp([tab, 'equilibrate = (none)']) ;
    end
    disp([tab, 'C = ', num2str(C)]) ;
    disp([tab, 'maxsteps = ', num2str(maxsteps)]) ;
    if hasTinitMethod
        disp([tab, 'Tinit = ', func2str(Tinit)]) ;
    else
        disp([tab, 'Tinit = ', num2str(Tinit)]) ;
    end
    disp([tab, 'r = ', num2str(r)]) ;
    if hasTfinalMethod
        disp([tab, 'Tfinal = ', func2str(Tfinal)]) ;
    else
        disp([tab, 'Tfinal = ', num2str(Tfinal)]) ;
    end
    disp([tab, 'f = ', num2str(f)]) ;
    disp([tab, 'maxtemps = ', num2str(maxtemps)]) ;
    disp([tab, 'v = ', num2str(v)]) ;
    disp([tab, 'bins = ', num2str(bins)]) ;
    disp([tab, 'e = ', num2str(e), newline]) ;
end
%
% Perform temperature initialization (temperature step 0).
%
if hasTinitMethod
    [T,W,Ew,Wbsf,Ebsf,Ea,Ev,steps] = feval(Tinit,r, walkers, newstate, X, cost, moveclass) ;
else
    [T,W,Ew,Wbsf,Ebsf,Ea,Ev,steps] = TinitT0(Tinit, walkers, newstate, X, cost, moveclass) ;
end
%
% Initialize counters, histories, etc.
% Note: Matlab matrix indicies run from 1 to whatever.
%   Consequently, temperature steps 0 to maxsteps are
%   internally indexed from 1 to maxsteps+1.
%
Eh = historyupdate([],Ev,0,Inf) ;
j = 1 ;
Tt(j) = Inf ;
Et(j) = Ea ;
Etarget(j) = Ea ;
ert(j) = 0 ;
Kt(j) = steps ;
Ebsft(j) = min(Ebsf) ;
%
if verbose
    disp(sprintf('%8s %10s %10s %10s %10s %10s %10s %10s %12s','t','T','<E>','Etarget', 'Estd','e','steps','Ebsf')) ;
    disp(sprintf('%8d %10.3g %10.3g %10.3g %10.3g %10.3g %10d %12.5g', ...
        round(j-1),Tt(j), Et(j), Etarget(j), std(Ew), ert(j), round(Kt(j)), Ebsft(j))) ;
end
%
% Iterature through the temperature steps
%
for i=1:maxtemps
    clear Ev ;
    %
    % Go take an equilibrium walk
    %
    [W,Ew,Wbsf,Ebsf,Ea,Estd,Ev,steps] = metropoliswalk( ...
        verbose, ...
        Ea, T, ...
        walkers, W, X, cost, moveclass, ...
        acceptrule, q, ...
        hasEquilibrate, equilibrate, C, maxsteps, ...
        Wbsf, Ebsf) ;
    %
    % Record what happenned
    %
    j = i+1 ;
    Tt(j) = T ;
    Et(j) = Ea ;
    Etarget(j) = Ea - (v*Estd) ;   % update the target on-the-fly
    Kt(j) = steps ;
    Ebsft(j) = min(Ebsf) ;
    Eh = historyupdate(Eh,Ev,i,T) ;  % update the history array
    %
    % Compute the density of states and relaxation time
    %
    if bins <= 0    % unless turned off
        M = [] ;
        rho = [] ;
        Ebin = [] ;
    else
        [M, e, rho, Ebin] = TM(Eh,bins) ;
    end
    ert(j) = e ;    % record the relaxation time
    %
    if verbose
        disp(sprintf('%8d %10.3g %10.3g %10.3g %10.3g %10.3g %10d %12.5g', ...
            round(j-1),Tt(j), Et(j), Etarget(j), Estd, ert(j), round(Kt(j)), Ebsft(j))) ;
    end
    %
    % Perform the temperature update.  Halt if some stopping criteria reached.
    %
    dEtgt = Etarget(j) - Etarget(j-1) ;
    T = feval(schedule,Ea,Estd,walkers,dEtgt,v,e,T,i,P) ;
    if hasTfinalMethod
        if feval(Tfinal,W,Ew,j,Tt,Et,Etarget,ert,Kt,Ebsft,f)
            if verbose
                disp(tab) ;
                disp([tab,'Stop criteria met for "',func2str(Tfinal),'"']) ;
            end
            break ;
        end
    elseif T < Tfinal
        if verbose
            disp(tab) ;
            disp([tab,'Tfinal (',num2str(Tfinal),') surpassed at T = ',num2str(T)]) ;
        end
        break ;
    end
end
if verbose
   disp(tab) ;
end
%
% Return data to caller.
%