initialize_sa Subroutine

   subroutine initialize_sa(me,n,lb,ub,c,&
                            maximize,eps,ns,nt,neps,maxevl,&
                            iprint,iseed1,iseed2,step_mode,vms,iunit,&
                            use_initial_guess,n_resets,&
                            cooling_schedule,cooling_param,&
                            optimal_f_specified,optimal_f,optimal_f_tol,&
                            distribution_mode,dist_std_dev,&
                            dist_scale,dist_shape,&
                            fcn,n_inputs_to_send,fcn_parallel_input,fcn_parallel_output,&
                            ireport,report)

      class(simulated_annealing_type),intent(inout) :: me
      integer, intent(in)                           :: n
      real(wp), dimension(n), intent(in)            :: lb
      real(wp), dimension(n), intent(in)            :: ub
      real(wp), dimension(n), intent(in),optional   :: c
      logical, intent(in),optional                  :: maximize
      real(wp), intent(in),optional                 :: eps
      integer, intent(in),optional                  :: ns
      integer, intent(in),optional                  :: nt
      integer, intent(in),optional                  :: neps
      integer, intent(in),optional                  :: maxevl
      integer, intent(in),optional                  :: iprint
      integer, intent(in),optional                  :: iseed1
      integer, intent(in),optional                  :: iseed2
      integer,intent(in),optional                   :: step_mode
      real(wp), intent(in), optional                :: vms
      integer, intent(in), optional                 :: iunit
      logical, intent(in), optional                 :: use_initial_guess
      integer, intent(in), optional                 :: n_resets
      integer, intent(in), optional                 :: cooling_schedule   !! temperature reduction schedule (1-5)
      real(wp), intent(in), optional                :: cooling_param      !! parameter for cooling schedules 3 and 5
      logical, intent(in), optional  :: optimal_f_specified !! if the optional `f` value is known,
                                                            !! it can be specified by `optimal_f`.
                                                            !! [Default is False]
      real(wp), intent(in), optional  :: optimal_f          !! if `optimal_f_specified=True` the solver will stop
                                                            !! if this value is achieved.
      real(wp), intent(in), optional  :: optimal_f_tol      !! absolute tolerance for the `optimal_f` check
      integer, dimension(:), intent(in), optional   :: distribution_mode  !! distribution for perturbations (per 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)
                                                                          !!
                                                                          !! Can be scalar (broadcast to all) or array(n)
      real(wp), dimension(:), intent(in), optional  :: dist_std_dev       !! std dev for normal/truncated_normal (per variable or scalar)
      real(wp), dimension(:), intent(in), optional  :: dist_scale         !! scale for cauchy/pareto (per variable or scalar)
      real(wp), dimension(:), intent(in), optional  :: dist_shape         !! shape for pareto (per variable or scalar)
      procedure(sa_func),optional                       :: fcn                  !! the user's function for scalar evaluation.
                                                                                !! if not present, the user must provide all of
                                                                                !! n_inputs_to_send, fcn_parallel_input and fcn_parallel_output
                                                                                !! for parallel evaluation.
      procedure(sa_func_parallel_inputs),optional       :: n_inputs_to_send     !! 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),optional  :: fcn_parallel_input   !! this function receives the input points for evaluation
      procedure(sa_func_parallel_output_func),optional  :: fcn_parallel_output  !! 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.
      integer, intent(in), optional      :: ireport    !! 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),optional :: report     !! if associated, this function is called to report intermediate
                                                       !! results to the user.

      if (present(fcn)) then
         ! serial function evaluation:
         me%fcn => fcn
         me%parallel_mode = .false.
      else if (present(n_inputs_to_send) .and. &
         present(fcn_parallel_input) .and. &
         present(fcn_parallel_output)) then
         ! parallel function evaluation:
         me%n_inputs_to_send    => n_inputs_to_send
         me%fcn_parallel_input  => fcn_parallel_input
         me%fcn_parallel_output => fcn_parallel_output
         me%parallel_mode = .true.
      else
         error stop 'Error: either fcn (serial mode) or all of n_inputs_to_send, '//&
                    'fcn_parallel_input and fcn_parallel_output (parallel mode) must be provided.'
      end if

      ! function evaluation reporting:
      if (present(report)) me%report => report
      if (present(ireport)) me%ireport = ireport

      me%n  = n
      me%lb = lb
      me%ub = ub

      ! check validity of bounds:
      if (any(me%lb > me%ub)) then
         error stop 'Error: at least one LB > UB.'
      end if

      allocate(me%c(n))
      if (present(c)) then
         me%c = c
      else
         me%c = 2.0_wp
      end if

      if (present(maximize)            ) me%maximize = maximize
      if (present(eps)                 ) me%eps = abs(eps)
      if (present(ns)                  ) me%ns = ns
      if (present(nt)                  ) me%nt = nt
      if (present(neps)                ) me%neps = neps
      if (present(maxevl)              ) me%maxevl = maxevl
      if (present(iprint)              ) me%iprint = iprint
      if (present(iseed1)              ) me%iseed1 = iseed1
      if (present(iseed2)              ) me%iseed2 = iseed2
      if (present(step_mode)           ) me%step_mode = step_mode
      if (present(vms)                 ) me%vms = vms
      if (present(iunit)               ) me%iunit = iunit
      if (present(use_initial_guess)   ) me%use_initial_guess = use_initial_guess
      if (present(n_resets)            ) me%n_resets = n_resets
      if (present(cooling_schedule)    ) me%cooling_schedule = cooling_schedule
      if (present(cooling_param)       ) me%cooling_param = cooling_param
      if (present(optimal_f_specified) ) me%optimal_f_specified = optimal_f_specified
      if (present(optimal_f)           ) me%optimal_f = optimal_f
      if (present(optimal_f_tol)       ) me%optimal_f_tol = abs(optimal_f_tol)

      ! allocate and set distribution parameters (per-variable):
      allocate(me%distribution_mode(n))
      allocate(me%dist_std_dev(n))
      allocate(me%dist_scale(n))
      allocate(me%dist_shape(n))

      ! set default values (uniform distribution with default parameters):
      me%distribution_mode = sa_mode_uniform
      me%dist_std_dev = 1.0_wp
      me%dist_scale = 1.0_wp
      me%dist_shape = 1.0_wp

      ! override with user-supplied values (broadcast if scalar):
      if (present(distribution_mode)) then
         if (size(distribution_mode) == 1) then
            me%distribution_mode = distribution_mode(1)  ! broadcast scalar
         else if (size(distribution_mode) == n) then
            me%distribution_mode = distribution_mode
         else
            error stop 'Error: distribution_mode must be scalar or size n.'
         end if
      end if
      if (present(dist_std_dev)) then
         if (size(dist_std_dev) == 1) then
            me%dist_std_dev = abs(dist_std_dev(1))
         else if (size(dist_std_dev) == n) then
            me%dist_std_dev = abs(dist_std_dev)
         else
            error stop 'Error: dist_std_dev must be scalar or size n.'
         end if
      end if
      if (present(dist_scale)) then
         if (size(dist_scale) == 1) then
            me%dist_scale = abs(dist_scale(1))
         else if (size(dist_scale) == n) then
            me%dist_scale = abs(dist_scale)
         else
            error stop 'Error: dist_scale must be scalar or size n.'
         end if
      end if
      if (present(dist_shape)) then
         if (size(dist_shape) == 1) then
            if (dist_shape(1) <= 0.0) then
               error stop 'Error: dist_shape must be strictly positive.'
            end if
            me%dist_shape = abs(dist_shape(1))
         else if (size(dist_shape) == n) then
            if (any(dist_shape <= 0.0)) then
               error stop 'Error: all dist_shape values must be strictly positive.'
            end if
            me%dist_shape = abs(dist_shape)
         else
            error stop 'Error: dist_shape must be scalar or size n.'
         end if
      end if

      ! validate distribution modes:
      if (any(me%distribution_mode < 1 .or. me%distribution_mode > n_distribution_modes)) then
         error stop 'Error: invalid distribution_mode.'
      end if

   end subroutine initialize_sa