module tspevol ! ! Evolutionary method for solving TSP (like Fogel, 1988) ! use tspglobals use tsputils use quicksort implicit none contains subroutine run_evolutionary integer, allocatable, dimension(:,:) :: parents, offspring real, allocatable, dimension(:) :: parents_length, offspring_length real, allocatable, dimension(:) :: pop_fitness, pop_indices real :: distance, new_distance real :: tem ! temperature integer :: tem_change_index, tem_trial_index integer :: i1, i2, i3 integer :: last_sample_trial integer :: gens integer :: n_compete integer :: n_offspring, b, i, j, uniq, diversity real :: rand print *, 'Using Evolutionary method.' n_offspring = birth_rate * n_parents + n_parents n_compete = n_offspring / 10 ! compete with 10% in tournament allocate(parents(0:n_parents-1, 0:n_cities-1)) allocate(offspring(0:n_offspring-1, 0:n_cities-1)) allocate(parents_length(0:n_parents-1), offspring_length(0:n_offspring-1)) allocate(pop_fitness(1:n_offspring), pop_indices(1:n_offspring)) ! initial tour if (init_type == INIT_NUMERICAL) then do i = 0, n_parents-1 call numerical_tour(parents(i,:)) end do else if (init_type == INIT_RANDOM) then do i = 0, n_parents-1 call random_tour(parents(i,:)) end do else if (init_type == INIT_GREEDY) then do i = 0, n_parents-1 call greedy_tour(parents(i,:)) end do end if ! calculate distances do i = 0, n_parents-1 parents_length(i) = tour_len(parents(i,:)) ! evaluate if (parents_length(i) < distance_min) then distance_min = parents_length(i) city_order_min = parents(i,:) end if end do ! other initialisation gens = 1 trials = 1 last_sample_trial = 1 tem = tem_init print *, 'trials ', 'distance_min ', 'tem ', 'gens ' print *, trials, distance_min, tem, gens 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 ! form a new generation !!do i = 0, n_parents-1 !! print *, 'parent', i, parents(i,:) !!end do loop_parents: do i = 0, n_parents-1 loop_birth: do b = 0, birth_rate-1 city_order = parents(i,:) distance = parents_length(i) ! mutate if (perturb_type == PERTURB_SWAP) then call random_number(rand) i1 = rand * n_cities call random_number(rand) i2 = rand * n_cities distance = swapped_distance(city_order, distance, i1, i2) call swap_operator(city_order, i1, i2) 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 distance = reversed_distance(city_order, distance, i1, i2) call reverse_operator(city_order, i1, i2) else ! shift segment distance = shifted_distance(city_order, distance, i1, i2, i3) call shift_operator(city_order, i1, i2, i3) end if end if offspring(i*birth_rate + b,:) = city_order offspring_length(i*birth_rate + b) = distance end do loop_birth ! include parents offspring(n_parents*birth_rate + i,:) = parents(i,:) offspring_length(n_parents*birth_rate + i) = parents_length(i) end do loop_parents ! calculate fitness ! pop_indices = 0; ! do i = 0, n_offspring-1 ! pop_indices(i+1) = real(i) ! print *, i, real(i), pop_indices(i+1) ! end do if (select_type == SELECT_RANK) then ! simple rank selection pop_indices = (/ (i,i=0,n_offspring-1) /) pop_fitness(1:n_offspring) = offspring_length(0:n_offspring-1) call sort2(n_offspring, pop_fitness, pop_indices) else if (select_type == SELECT_BOLTZ_10) then ! Boltzmann selection, 10% competition pop_indices = (/ (i,i=0,n_offspring-1) /) pop_fitness(1:n_offspring-1) = 0 do i = 0, n_offspring-1 ! compete against 10% do j = 0, n_compete call random_number(rand) i2 = rand * n_offspring if (metrop_boltz(offspring_length(i) - offspring_length(i2), tem)) then pop_fitness(i+1) = pop_fitness(i+1) + 1 else ! not necessary, but allows 'upsets' pop_fitness(i2+1) = pop_fitness(i2+1) + 1 end if end do end do call sort2(n_offspring, pop_fitness, pop_indices) ! reverse so that highest fitness is at start of array !pop_fitness = pop_fitness(n_cities:1:-1) !pop_indices = pop_indices(n_cities:1:-1) end if !print *, 'fitness: ', pop_fitness !print *, 'indices: ', pop_indices ! select survivors do i = 0, n_parents-1 j = pop_indices(i+1) if (select_type == SELECT_BOLTZ_10) then ! fitness sorted ascending; select from end j = pop_indices(n_offspring - i) end if parents(i,:) = offspring(j,:) parents_length(i) = offspring_length(j) ! evaluate if (parents_length(i) < distance_min) then distance_min = parents_length(i) city_order_min = parents(i,:) if (show_chgs) then if (trials - last_sample_trial >= sampling_threshold) then print *, trials, distance_min, tem, gens last_sample_trial = trials end if end if end if end do trials = trials + n_parents*birth_rate gens = gens + 1 ! recalculate true distance occasionally to avoid rounding errors if (mod(gens, recalc_trials) == 0) then do i = 0, n_parents-1 parents_length(i) = tour_len(parents(i,:)) end do end if ! print status out_samples times if (mod(gens, out_sample_trials) == 0) then print *, trials, distance_min, tem, gens end if end do loop_at_tem end do loop_schedule ! recalculate true distance do i = 0, n_parents-1 parents_length(i) = tour_len(parents(i,:)) end do ! calculate diversity pop_fitness(1:n_parents) = parents_length(0:n_parents-1) call sort2(n_parents, pop_fitness, pop_indices) uniq = 1 do i = 2, n_parents if (pop_fitness(i) /= pop_fitness(i-1)) then uniq = uniq + 1 end if end do diversity = 100.0 * real(uniq) / real(n_parents) print *, '' print *, 'Diversity: ', diversity, '%' end subroutine run_evolutionary end module tspevol