anneal.f90.html

!
! 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