module tspanneal

  !
  ! Simulated Annealing method for solving TSP (Boltzmann selection)
  !

  use tspglobals
  use tsputils

  implicit none

  contains

    subroutine run_annealing
      real :: distance, new_distance
      real :: tem ! temperature
      integer :: tem_change_index, tem_trial_index
      integer :: tour_changes
      integer :: i1, i2, i3
      integer :: last_sample_trial
      real :: rand

      print *, 'Using Simulated Annealing method.'

      ! initial tour
      if (init_type == INIT_NUMERICAL) then
         call numerical_tour(city_order)
      
      else if (init_type == INIT_RANDOM) then
         call random_tour(city_order)
         
      else if (init_type == INIT_GREEDY) then
         call greedy_tour(city_order)
      end if

      ! other initialisation
      trials = 1
      last_sample_trial = 1
      distance = tour_len(city_order)
      city_order_min = city_order
      distance_min = distance
      tem = tem_init
      tour_changes = 0

      print *, 'trials ', 'distance_min ', 'tem ', 'tour_changes '
      print *, trials, distance_min, tem, tour_changes

      loop_schedule: do tem_change_index = 1, sched_changes

         if (sched_type == SCHED_LINEAR) then
            tem = tem_init * (1 - real(tem_change_index) / real(sched_changes))

         else if (sched_type == SCHED_QUADRATIC) then
            tem = tem_init - sched_quad_const * tem_change_index**2

         else if (sched_type == SCHED_STEP) then
            !if (tem_change_index == sched_changes/2) then
            !   tem = 0
            !end if
            if (tem_change_index > sched_changes/2) then
               tem = 0
            else
               tem = tem_init * (1 - real(tem_change_index * 2) / real(sched_changes))
            end if
         end if

         loop_at_tem: do tem_trial_index = 1, trials_at_tem

            ! perturb tour
            if (perturb_type == PERTURB_SWAP) then
               call random_number(rand)
               i1 = rand * n_cities
               call random_number(rand)
               i2 = rand * n_cities
               new_distance = swapped_distance(city_order, distance, i1, i2)
               ! determine whether to accept perturbation
               if (metrop_boltz(new_distance - distance, tem)) then
                  call swap_operator(city_order, i1, i2)
                  distance = new_distance
                  tour_changes = tour_changes + 1
               end if

            else if (perturb_type == PERTURB_REVSHIFT) then
               call random_number(rand)
               i1 = rand * n_cities
               call random_number(rand)
               i2 = rand * n_cities
               ! reverse or shift?
               call random_number(rand)
               if (rand < 0.5) then
                  ! reverse segment
                  new_distance = reversed_distance(city_order, distance, i1, i2)
                  ! determine whether to accept perturbation
                  if (metrop_boltz(new_distance - distance, tem)) then
                     call reverse_operator(city_order, i1, i2)
                     distance = new_distance
                     tour_changes = tour_changes + 1
                  end if
               else
                  ! shift segment
                  new_distance = shifted_distance(city_order, distance, i1, i2, i3)
                  ! determine whether to accept perturbation
                  if (metrop_boltz(new_distance - distance, tem)) then
                     call shift_operator(city_order, i1, i2, i3)
                     distance = new_distance
                     tour_changes = tour_changes + 1
                  end if
               end if
            end if
            
            trials = trials + 1

            ! recalculate true distance occasionally to avoid rounding errors
            if (mod(trials, recalc_trials) == 0) then
               distance = tour_len(city_order)
            end if

            ! evaluate
            if (distance < distance_min) then
               distance_min = distance
               city_order_min = city_order
               if (show_chgs) then
                  if (trials - last_sample_trial >= sampling_threshold) then
                     print *, trials, distance_min, tem, tour_changes
                     last_sample_trial = trials
                  end if
               end if
            end if

            ! print status out_samples times
            if (mod(trials, out_sample_trials) == 0) then
               print *, trials, distance_min, tem, tour_changes
            end if
            
         end do loop_at_tem
      end do loop_schedule

    end subroutine run_annealing

end module tspanneal