!
! simulated annealing to solve problem 4 from pt II course
!
! copyright (c) 1998 Chris Lightfoot
! cwrl2@hermes.cam.ac.uk
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! GLOBAL TYPES
!
MODULE types
TYPE node_pos
REAL:: x, y;
END TYPE node_pos;
END MODULE types
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! SETUP SUBROUTINES
!
! read data from file into an array of type node_pos; the input file is
! assumed to be formatted as /([0-9.]{7}) +([0-9.]{7})/ or similar
SUBROUTINE read_data(filename, num_nodes, node)
USE types
IMPLICIT NONE
CHARACTER (LEN=*), INTENT(in):: filename
INTEGER, INTENT(in):: num_nodes
TYPE(node_pos), DIMENSION(num_nodes):: node
INTEGER:: i
! uncomment this to print out the filename read
!
! PRINT 4211, filename
!4211 FORMAT("opening file ", A);
OPEN(unit=42, file=filename, action="read", err=4266)
DO i = 1, num_nodes
READ(unit=42, fmt=4200) node(i)%x, node(i)%y
4200 FORMAT(2f7.3)
! uncomment this to print out the node coordinates
!
! PRINT 4201, i, node(i)%x, node(i)%y
!4201 FORMAT(i4, " [", f5.3, "] [", f5.3, "]");
END DO
CLOSE(42);
RETURN
4266 CALL PERROR("read_data") ! perror is defined in stupid.a
STOP "an error occurred: bailing out"
END SUBROUTINE read_data
! calculate internodal distances, placing them in a large matrix for
! later reference
SUBROUTINE calculate_distances(num_nodes, nodes, distances)
USE types
IMPLICIT NONE
INTEGER, INTENT(in):: num_nodes
TYPE(node_pos), DIMENSION(num_nodes), INTENT(IN):: nodes
REAL, DIMENSION(num_nodes, num_nodes), INTENT(OUT):: distances
INTEGER:: i, j
REAL:: s
DO i = 1, num_nodes
distances(i, i) = 0;
! the matrix of distances is symmetric about its leading diagonal; we
! exploit this in order to evaluate only about 1/2 * num_nodes ** 2
! square roots
DO j = 1, i - 1
s = SQRT((nodes(i)%x - nodes(j)%x)**2 + (nodes(i)%y - nodes(j)%y)**2);
distances(i, j) = s;
distances(j, i) = s;
END DO
END DO
END SUBROUTINE calculate_distances
! retrieve the parameters for the simulation, which are (in order) the
! iterations_limit, the initial alpha and the cooling factor (as a
! fraction of alpha)
SUBROUTINE get_parameters(iterations_limit, alpha, cooling, swap_prob)
INTEGER, INTENT(OUT):: iterations_limit
REAL, INTENT(OUT):: alpha, cooling, swap_prob
! uncomment this to prompt for parameters
!
! WRITE(*,*) "enter iterations_limit, alpha, cooling, swap_prob"
READ(*, 5000) iterations_limit
5000 FORMAT(i8)
READ(*, *) alpha
READ(*, *) cooling
READ(*, *) swap_prob
! uncomment this to print out the parameters given
!
! WRITE(*, 5002) iterations_limit, alpha, cooling, swap_prob
!5002 FORMAT("iterations_limit = ", i8, "; alpha = ", f9.7, "; cooling = ", f9.7, "; swap_prob = ", f9.7)
END SUBROUTINE get_parameters
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! NUMERICAL ROUTINES
!
! given an integer array representing a configuration, calculate the
! path length of the configuration using a matrix of internodal
! distances calculated above
REAL FUNCTION find_config_length(num_nodes, config, distances)
IMPLICIT NONE
INTEGER, INTENT(IN):: num_nodes
INTEGER, DIMENSION(num_nodes), INTENT(IN):: config
REAL, DIMENSION(num_nodes, num_nodes), INTENT(in):: distances
REAL:: s
INTEGER:: i
! initialise the path length to 0
s = 0
DO i = 1, (num_nodes - 1)
s = s + distances(config(i), config(i + 1))
END DO
find_config_length = s
END FUNCTION find_config_length
! this is not really a numerical routine, it's just here to hide the
! evil f77 NAG libraries from my nice post-stone-age code
INTEGER FUNCTION random_integer(max)
IMPLICIT NONE
INTEGER, INTENT(IN):: max
INTEGER:: G05DYF ! you must be joking
! G05DYF (how /do/ they come up with these names?) returns an integer
! formed from
!
! (min value) + [(n - m + 1)y] ([] denote the integer part)
!
! so, since we are interested in numbers in [1, max] as array indices,
! we pass 1 and max as the parameters
random_integer = G05DYF(1, max)
END FUNCTION random_integer
! similarly, to wrap NAG's G05CAF (real uniform deviates over [0,1])
REAL FUNCTION random_real()
IMPLICIT NONE
REAL:: x
REAL:: G05CAF ! bletch
! G05CAF is the "generic" uniform deviates random number generator
random_real = G05CAF(x) ! x does nothing but is required by the Standard
END FUNCTION random_real
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! UTILITY SUBROUTINES
!
! print out a path configuration and its length
SUBROUTINE print_config(num_nodes, config, len)
IMPLICIT NONE
INTEGER, INTENT(IN):: num_nodes
INTEGER, DIMENSION(num_nodes), INTENT(IN):: config
REAL, INTENT(IN):: len
! REAL:: find_config_length
INTEGER:: i
WRITE(unit=*, fmt=4201, advance='no') len
4201 FORMAT(1X, F7.3)
DO i = 1, num_nodes
WRITE(unit=*, fmt=4200, advance='no') config(i)
4200 FORMAT(1X, I3)
END DO
PRINT *
END SUBROUTINE print_config
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! MAIN PROGRAM
!
PROGRAM anneal
USE types
IMPLICIT NONE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! FUNCTION DECLARATIONS
!
INTEGER:: random_integer
REAL:: random_real
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! GLOBAL DATA
!
! number of nodes
INTEGER, PARAMETER:: num_nodes = 40
! positions of nodes
TYPE(node_pos), DIMENSION(num_nodes):: node_list
! internodal distances
REAL, DIMENSION(num_nodes, num_nodes):: distance_matrix
! configuration of loop
INTEGER, DIMENSION(num_nodes):: config, config_new, config_best
REAL:: length, length_new, length_best
! random index variables
INTEGER:: i
! functions
REAL:: find_config_length
!
! parameters to do with the annealing
!
! how long we go on for if there are no configuration changes
INTEGER:: iterations_limit
! where we've got to
INTEGER:: iterations = 0
! the thermodynamic temperature (alpha == something like 1 / kT)
REAL:: alpha
! cooling factor per iteration
REAL:: cooling
! the probability that we swap two nodes, rather than reversing the order of
! a section of the path
REAL:: swap_prob
! stuff to do with reconfigurations
INTEGER:: a, b, c
! initialise the random number generator
CALL G05CCF ! bletch
! retrieve information about how to run the simulation
CALL get_parameters(iterations_limit, alpha, cooling, swap_prob);
! read in the data, and process it to generate a matrix of distances
CALL read_data("station.data", num_nodes, node_list);
CALL calculate_distances(num_nodes, node_list, distance_matrix);
! establish an initial configuration
!
! it doesn't matter what the configuration is, since early in the annealing
! process, total randomisation of the system will take place. hence, we
! choose the simplest possible one
DO i = 1, num_nodes
config(i) = i
END DO
length = find_config_length(num_nodes, config, distance_matrix)
length_best = length
! uncomment this to print out the initial configuration
!
!CALL print_config(num_nodes, config, length)
! this is the main loop which runs the annealing algorithm
DO WHILE (iterations < iterations_limit)
! try a new configuration
config_new = config;
! now we try a reconfiguration
a = random_integer(num_nodes)
b = random_integer(num_nodes)
DO WHILE (b == a)
b = random_integer(num_nodes)
END DO
! select a reconfiguration type
IF (random_real() < swap_prob) THEN
! this is a reconfiguration based on swapping two nodes
c = config_new(a)
config_new(a) = config_new(b)
config_new(b) = c
ELSE
! this is a reconfiguration based on inverting the order of a section of
! the path
IF (b < a) THEN
c = a
a = b
b = c
END IF
DO c = a, b
config_new(c) = config(a + b - c)
END DO
END IF
! test the reconfiguration
length_new = find_config_length(num_nodes, config_new, distance_matrix)
! if the new configuration is shorter than the previous one, accept it;
! otherwise, accept it with probability exp((length - length_new) / kT)
IF (length_new < length) THEN
config = config_new
length = length_new
ELSE IF (random_real() < EXP((length - length_new) * alpha)) THEN
config = config_new
length = length_new
END IF
! now, if the new configuration is better than the best already-established
! configuration, we keep it and print it out
IF (length < length_best) THEN
config_best = config
length_best = length
iterations = 0
! uncomment this to print out the configuration every time a new better one
! is found
!
! CALL print_config(num_nodes, config, length)
END IF
! decrease the temperature a little
alpha = alpha + alpha * cooling
! increment iterations
iterations = iterations + 1
END DO
CALL print_config(num_nodes, config, length)
END PROGRAM anneal
Copyright (c) 1999 Chris Lightfoot. All rights reserved.