function [W,Ew,Wbsf,Ebsf,Tt,Et,Etarget,ert,Kt,Ebsft,Eh,M,rho,Ebin] = anneal( ...
    verbose, ...
    newstate, X, ...
    cost, moveclass, ...
    walkers, ...
    acceptrule, q, ...
    schedule, P, ...
    equilibrate, C, maxsteps, ...
    Tinit, r, ...
    Tfinal, f, maxtemps, ...
    v, bins, e)
% MAIN DRIVER and HELP file supplied with SA Tools.
% Copyright (c) 2002, by Richard Frost and Frost Concepts.
% See http://www.frostconcepts.com/software for information on SA Tools.
% Get the book:  http://www.frostconcepts.com/books/ebsa/
%
% [W,Ew,Wbsf,Ebsf,Tt,Et,Etarget,ert,Kt,Ebsft,Eh,M,rho,Ebin] = anneal( ...
%     verbose, ...
%     newstate, X, ...
%     cost, moveclass, ...
%     walkers, ...
%     acceptrule, q, ...
%     schedule, P, ...
%     equilibrate, C, maxsteps, ...
%     Tinit, r, ...
%     Tfinal, f, maxtemps, ...
%     v, bins, e)
%
%   verbose = prints status information when true (1).
%   newstate = (handle to) user-defined method
%           W0 = newstate(X)    where
%               X = user-defined problem domain or other data,
%                       behaviorally static.
%               W0 = an initial user-defined state.
%           Book chapter 2.
%   X = user-defined problem domain or other data, behaviorally static.
%           Book chapter 2.
%   cost = (handle to) user-defined objective method (function)
%           Ew = cost(X,W)    where
%               X = user-defined problem domain or other data.
%               W = a user-defined state from 'newstate' or 'moveclass'.
%               Ew = energy corresponding to W
%           Book chapter 9.
%   moveclass = (handle to) user-defined method,
%           W = moveclass(X,W,Ea,T)    where
%               X = user-defined problem domain or other data.
%               W = a user-defined state from 'newstate' or 'moveclass'.
%               Ea = average energy at current temperature.
%               T = current temperature
%           Book chapters 2.2 and 10.2.
%   walkers = number of walkers.  Must be positive integer.
%               walkers = 1 implies barebones annealing
%               walkers > 4 suggested for ensemble methods
%           Book chapters 4 and 7.
%   acceptrule = (handle to) SA Tools or user-defined method
%           a = acceptrule(dE,T,q)    where
%               dE = the difference in cost between a trial state and
%                       the current state: dE = Wtrial - W
%               T = the current temperature
%               q = any data required by the acceptrule
%               a = 0 if trial is rejected, otherwise 1.
%           SA Tools supplied methods are:
%               metropolis
%               szu
%               tsallis
%               threshold
%               franz
%           Book chapter 11.    
%   q = any data required by the acceptrule.
%           Book chapter 11.
%   schedule = (handle to) SA Tools or user-defined temperature update
%           nextT = schedule(Ea,Estd,walkers,dEtgt,v,e,T,t,P)    where
%               Ea = average energy at current temperature.
%               Estd = standard deviation of energies
%               dEtgt = difference between present and previous target mean energy
%               walkers = number of walkers.  Must be positive integer.
%               T = current temperature
%               i = # of current temperature
%                   (i.e., 1st temperature is 1, 2nd is 2, etc.)
%               P = any data required by schedule
%               nextT = next temperature
%           SA Tools supplied methods are:
%               geman
%               geometric
%               hartley
%               berkeley
%               thermospeedHC
%               thermospeedR
%               retrospect
%           Book chapter 13.
%   P = any data required by schedule.
%   equilibrate = (handle to) SA Tools method, or user-defined method,
%           or a non-function_handle type (e.g., 0).  If a function handle is 
%           supplied, then the temperature will not change (i.e., schedule will
%           not be called) until equilibrate returns false (0).  Otherwise, the
%           moveclass will be executed maxsteps times between each
%           temperature change.  Method signature:
%               b = equilibrate(Ea0,Ea,Ew,walkers,T,step,maxsteps,C)    where
%                   Ea0 = average energy at the beginning of the metropolis walk
%                   Ea = current average energy
%                   Ew = current energies corresponding to W (size walkers)
%                   walkers = the number of walkers in the simulation
%                   T = the current temperature
%                   step = the current number of steps taken in the walk
%                   maxsteps = an upper limit on the number of steps in the walk
%                   C = any behaviorally constant data required by the method
%                   b = 0 if the temperature may change, otherwise 1.
%               SA Tools supplied methods are:
%                   hoffmann    (wait-for-a-fluctuation)
%           Book chapter 13.
%   C = any data required by equilibrate.
%   maxsteps = maximum number of times to attempt equilibration (call moveclass at fixed T).
%           Book chapter 13.
%   Tinit = initial temperature (Inf ok) -- or (handle to) SA Tools method,
%           or user-defined method.  If method handle is not present, then the
%           initial temperature will be T0 = Tinit.  Otherwise, the method will
%           calculate T0.  All moves made during this method must be accepted.
%           Method signature:
%               [T0,W,Ew,Ev,steps] = Tinit(r, walkers, newstate, X, cost, moveclass)
%                   INPUTS:
%                       r = behaviorially constant data required by Tinit (if any)
%                       walkers, newstate, X, cost, moveclass: defined above
%                   OUTPUTS:
%                       T0 = initial temperature
%                       steps = # of steps taken by each walker during Tinit
%                       Ev = energy (cost) history at T (infinite for Tinit)
%                           i = arbitrary index
%                           Ev(i,1) = step #
%                           Ev(i,2) = walker #
%                           Ev(i,3) = an energy visited during T
%                           Ev(i,4) = energy attempted from Ev(i,1:3) during T
%                       W,Ew: defined below
%               SA Tools supplied methods are:
%                   TinitT0
%                   TinitAccept
%                   TinitWhite
%           Book section 13.1.
%   r = behaviorially constant data required by Tinit (if any)
%           Book section 13.1.
%   Tfinal = final temperature (-Inf ok) or (handle to) SA Tools method,
%           or user-defined method.  If method handle is not present, then the
%           simulation will end when T drops below the value of Tfinal.
%           Otherwise, the method will calculate a logical value which when
%           true (equal to 1) will stop the simulation.
%           Method signature:
%               b = Tfinal(W,Ew,t,Tt,Et,Etarget,ert,Kt,Ebsft,f)
%                   INPUTS:
%                       W = cell array of current states (size walkers)
%                       Ew = energies associated with W
%                       t = current temperature step index; i.e., current T = Tt(t).
%                       Tt = temperature history of simulation (so far)
%                       Et = mean energy history
%                       Etarget = target mean energy history
%                       ert = relaxation time history
%                       Kt = equilibrium step history
%                       Ebsft = Ebsf history
%                       f = behaviorally constant data required by Tfinal method
%                   OUTPUT:
%                       b = true (equal to 1) when final temperature iteration has been reached
%               SA Tools supplied methods are:
%                   TfinalNstep
%           Book section 13.1.
%   f = behaviorally constant data required by Tfinal method
%           Book section 13.1.
%   maxtemps = maximum number of temperature iterations.
%           Book chapter 13.
%   v = thermodynamic speed.
%               Effects thermospeed schedules.
%               Typically 0 < v < 1.  0 ok for non-thermospeed schedules.
%           Book chapter 13.
%   bins = # of bins to use in estimation of M, e, rho, and Ebin each temperature step.
%               If bins <= 0, then M, e, rho, and Ebin will not be calculated
%                   and the user-supplied constant value of e will be used each step.
%               If bins > 0, then the supplied value of e will be ignored and
%                   the TM method will be called to calculate M, e, rho, and Ebin.
%           Book section 12.2.1.
%   e = estimate of relaxation time.  See bins, above.
%           Book section 12.2.1.
%
%   RETURN VALUES:
%   W = cell array of final state(s) of size 'walkers'
%   Ew = array of final energies corresponding to W
%   Wbsf = array of best-so-far states of size 'walkers'
%   Ebsf = array of best-so-far energies
%   Tt(i) = temperature at temperature step i-1
%       NOTE: matlab does not permit indicies less than 1,
%               so step 0 is at i=1, etc.
%   Et(i) = average energy at Tt(i)
%   Etarget(i) = target mean energy at Tt(i), calculated with v.
%   ert(i) = estimated relaxation time at Tt(i), calculated with bins.
%   Kt(i) = number of equilibration steps taken at Tt(i)
%   Ebsft(i) = best-so-far energy at Tt(i)
%   Eh = energy and temperature history
%          i = 1, 1+(steps*walkers), etc.
%          Eh(i,1) = index t of temperature step
%          Eh(i,2) = T corresponding to t
%          Eh(i,3) = equilibrium step #j at T
%          Eh(i,4) = walker #k
%          Eh(i,5) = energy E visited by walker k at step j during T
%          Eh(i,6) = energy E' attempted from E by walker k at step j during T
%   M = final Transition Matrix (see book section 12.2.1).
%   rho = final estimate of equilibrium density of states
%   Ebin = energy bin centroids, min, and max
%          Ebin(1,:) are bin centroids.  Ebin(1,b) is the centroid for rho(b).
%          Ebin(2,:) are bin lower bounds
%          Ebin(3,:) are bin upper bounds
%
% Example uses of this driver can be found in the examples/ directory.
%
%   e.g.,
%
%   rand('state',sum(100*clock)) ;
%   [W,Ew,Wbsf,Ebsf,Tt,Et,Etarget,ert,Kt,Ebsft,Eh,M,rho,Ebin] = anneal( ...
%       1, ...
%       12, ...