simulated_annealing_type Derived Type

   type,public :: simulated_annealing_type

      private

      integer :: n = 0              !! number of variables in the function to be optimized.
      logical :: maximize = .false. !! denotes whether the function should be maximized or minimized.
                                    !! a true value denotes maximization while a false value denotes
                                    !! minimization.  intermediate output (see iprint) takes this into
                                    !! account.
      real(wp) :: eps = 1.0e-9_wp   !! error tolerance for termination. if the final function
                                    !! values from the last neps temperatures differ from the
                                    !! corresponding value at the current temperature by less than
                                    !! eps and the final function value at the current temperature
                                    !! differs from the current optimal function value by less than
                                    !! eps, execution terminates and ier = 0 is returned.
      integer  :: ns = 20           !! number of cycles.  after ns function evaluations, each element of
                                    !! vm is adjusted according to the input `step_mode`
                                    !! the suggested value is 20.
      integer  :: nt = 100          !! number of iterations before temperature reduction. after
                                    !! nt*ns function evaluations, temperature (t) is changed
                                    !! by the factor rt.  value suggested by corana et al. is
                                    !! max(100, 5*n).  see goffe et al. for further advice.
      integer  :: neps = 4          !! number of final function values used to decide upon
                                    !! termination.  see eps.  suggested value is 4.
      integer  :: maxevl = 10000    !! the maximum number of function evaluations.  if it is
                                    !! exceeded, ier = 1.

      logical :: use_initial_guess = .true. !! if false, the initial guess is ignored and a
                                            !! random point in the bounds is used for the first
                                            !! function evaluation

      integer :: n_resets = 2 !! number of times to run the main loop (must be >=1)

      real(wp), dimension(:),allocatable :: lb   !! the lower bound for the allowable solution variables.
      real(wp), dimension(:),allocatable :: ub   !! the upper bound for the allowable solution variables.
                                                 !! if the algorithm chooses x(i) < lb(i) or x(i) > ub(i),
                                                 !! i = 1, n, a point is from inside is randomly selected. this
                                                 !! this focuses the algorithm on the region inside ub and lb.
                                                 !! unless the user wishes to concentrate the search to a particular
                                                 !! region, ub and lb should be set to very large positive
                                                 !! and negative values, respectively.  note that the starting
                                                 !! vector x should be inside this region.  also note that lb and
                                                 !! ub are fixed in position, while vm is centered on the last
                                                 !! accepted trial set of variables that optimizes the function.
      real(wp), dimension(:),allocatable :: c  !! vector that controls the step length adjustment.  the suggested
                                               !! value for all elements is 2.0.
      integer :: iprint = 1     !! controls printing inside [[sa]]:
                                !!
                                !!  * 0 - nothing printed.
                                !!  * 1 - function value for the starting value and
                                !!    summary results before each temperature
                                !!    reduction. this includes the optimal
                                !!    function value found so far, the total
                                !!    number of moves (broken up into uphill,
                                !!    downhill, accepted and rejected), the
                                !!    number of out of bounds trials, the
                                !!    number of new optima found at this
                                !!    temperature, the current optimal x and
                                !!    the step length vm. note that there are
                                !!    n*ns*nt function evalutations before each
                                !!    temperature reduction. finally, notice is
                                !!    is also given upon achieveing the termination
                                !!    criteria.
                                !!  * 2 - each new step length (vm), the current optimal
                                !!    x (xopt) and the current trial x (x). this
                                !!    gives the user some idea about how far x
                                !!    strays from xopt as well as how vm is adapting
                                !!    to the function.
                                !!  * 3 - each function evaluation, its acceptance or
                                !!    rejection and new optima. for many problems,
                                !!    this option will likely require a small tree
                                !!    if hard copy is used. this option is best
                                !!    used to learn about the algorithm. a small
                                !!    value for maxevl is thus recommended when
                                !!    using iprint = 3.
                                !!
                                !! suggested value: 1
                                !! note: for a given value of iprint, the lower valued
                                !! options (other than 0) are utilized.
      integer :: iseed1 = 1234    !! the first seed for the random number generator.
      integer :: iseed2 = 5678    !! the second seed for the random number generator.
                                  !! different values for iseed1 and iseed2 will lead
                                  !! to an entirely different sequence of trial points
                                  !! and decisions on downhill moves (when maximizing).
                                  !! see goffe et al. on how this can be used to test
                                  !! the results of [[sa]].
      integer   :: step_mode = 1  !! how to vary `vm` after `ns` cycles.
                                  !!
                                  !! * 1 : original method: adjust vm so that approximately
                                  !!   half of all evaluations are accepted
                                  !! * 2 : keep vm constant
                                  !! * 3 : adjust by a constant `vms` factor
      real(wp)  :: vms = 0.1_wp         !! for `step_mode=3`, the factor to adjust `vm`
      integer   :: iunit = output_unit  !! unit number for prints.

      ! cooling schedule options:
      integer   :: cooling_schedule = 1 !! temperature reduction schedule:
                                        !!
                                        !! * 1 : geometric (default): T(k+1) = rt * T(k)
                                        !! * 2 : fast annealing (Cauchy): T(k) = T0 / (1 + k)
                                        !! * 3 : huang: T(k) = T0 / (1 + cooling_param * k)^n
                                        !! * 4 : boltzmann: T(k) = T0 / log(1 + k + exp(1))
                                        !! * 5 : logarithmic: T(k) = T0 / (1 + cooling_param * log(1 + k))
      real(wp)  :: cooling_param = 1.0_wp !! parameter for cooling schedules 3 and 5. suggested value: 1.0

      ! option to specify the known answer, so the solver will stop if it finds it:
      logical :: optimal_f_specified = .false. !! if the optional `f` value is known,
                                               !! it can be specified by `optimal_f`.
      real(wp) :: optimal_f = 0.0_wp !! if `optimal_f_specified=True` the solver will stop
                                     !! if this value is achieved.
      real(wp) :: optimal_f_tol = 0.0_wp !! absolute tolerance for the `optimal_f` check

      ! distribution selection for perturbations (per-variable):
      integer, dimension(:), allocatable :: distribution_mode !! distribution to use for perturbations for each variable:
                                                              !!
                                                              !! * `sa_mode_uniform` : uniform (default)
                                                              !! * `sa_mode_normal` : normal (Gaussian)
                                                              !! * `sa_mode_cauchy` : cauchy
                                                              !! * `sa_mode_triangular` : triangular
                                                              !! * `sa_mode_bipareto` : bipareto (two-sided pareto)

      ! distribution parameters (per-variable, used depending on distribution_mode):
      real(wp), dimension(:), allocatable :: dist_std_dev   !! standard deviation for normal/truncated_normal
      real(wp), dimension(:), allocatable :: dist_scale     !! scale parameter for cauchy/bipareto distributions
      real(wp), dimension(:), allocatable :: dist_shape     !! shape parameter for bipareto distribution

      ! serial function evaluation:
      procedure(sa_func),pointer :: fcn => null()  !! the user's function

      ! function to report intermediate results to the user:
      integer :: ireport = 0 !! how often to report intermediate results to the user via `report` function:
                            !!
                            !! * 0 : no intermediate reports
                            !! * 1 : report each function evaluation
                            !! * 2 : report after each new optimal value is found
                            !! * 3 : report each function evaluation and each new optimal value found
      procedure(sa_report_func),pointer :: report => null() !! if associated, this function is called to report intermediate results to the user.

      ! parallel function evaluation (optional):
      logical :: parallel_mode = .false. !! if true, the user wants to evaluate multiple points in parallel (e.g. on a GPU or with OpenMP). if false, the user will evaluate one point at a time. if true, the user must provide all of n_inputs_to_send, fcn_parallel_input and fcn_parallel_output.
      procedure(sa_func_parallel_inputs),pointer :: n_inputs_to_send => null()  !! if the user wants to evaluate multiple points in parallel, this function should return the number of input points that will be sent for evaluation at a time.
      procedure(sa_func_parallel_inputs_func),pointer :: fcn_parallel_input => null()  !! this function receives the input points for evaluation
      procedure(sa_func_parallel_output_func),pointer :: fcn_parallel_output => null()  !! this function receives the output points from the parallel evaluation. The `x` here will be one of the ones sent via `sa_func_parallel_inputs_func`, The algorithm will only process one at the time.

      ! termination safeguards (to prevent freezing/hanging):
      real(wp) :: vm_min = 1.0e-12_wp  !! minimum allowable step size. if all vm values
                                       !! fall below this, terminate with ier=6

      contains

      private

      procedure,public :: initialize => initialize_sa
      procedure,public :: optimize   => sa
      procedure,public :: destroy    => destroy_sa

      procedure :: func
      procedure :: perturb_and_evaluate
      procedure,public :: perturb_variable  !! this is public so we can use it in the tests

   end type simulated_annealing_type