; XEmacs: This file contains -*-Scheme-*- source code.

;;; John David Stone
;;; Department of Mathematics and Computer Science
;;; Grinnell College
;;; stone@cs.grinnell.edu

;;; created December 5, 1999
;;; last revised December 6, 2000

;;; Simulated annealing

;;; When no efficient algorithm for solving an optimization problem is
;;; available, simulated annealing can often be used to find a near-optimal
;;; solution.  (However, it is not guaranteed to work; the justification
;;; for using it is entirely pragmatic.)

;;; The metaphor underlying this approach is the physical process of making
;;; a metal less brittle by annealing it -- heating it up and then cooling
;;; it very slowly, waiting at each cooling step for the atoms in the metal
;;; to find optimal positions and orientations, so as to reach the state of
;;; minimum energy for the system.

;;; The SIMULATED-ANNEALING procedure encapsulates the control structure
;;; for such a simulation.  Initially, one possible solution is constructed
;;; arbitrarily, perhaps randomly.  It is expected that this initial
;;; solution is far from optimal, so the solution is then repeatedly
;;; modified.  The successive modifications are also chosen randomly, but
;;; not all of the modifications are actually made.  Instead, the effect of
;;; each proposed modification on the quantity to be optimized is
;;; evaluated.  The proposed modification is actually made whenever it
;;; increases the value of the solution, but also sometimes when it
;;; decreases it by a sufficiently small amount.

;;; In early stages, even modifications that considerably decrease the
;;; value of the solution are sometimes approved, just as in early stages
;;; of annealing the atoms are relatively free to shift their positions and
;;; orientations.  But as the simulated temperature is lowered, the
;;; probability of a large decrease in the value of a solution is also
;;; gradually lowered and the pieces of the solution are therefore
;;; increasingly locked into place.

;;; The SIMULATED-ANNEALING procedure takes five arguments:

;;; INITIAL-SOLUTION is an arbitrarily constructed initial solution.

;;; The TRIAL-LIMIT parameter puts an upper bound on the number of
;;; modifications that will be proposed at any single temperature level.
;;; When this limit is reached, the temperature is lowered and the next
;;; phase of annealing begins.

;;; The CHANGE-LIMIT parameter puts an upper bound on the number of
;;; modifications that will actually be made at any single temperature
;;; level.  Reaching this limit also causes a transition to the next phase
;;; of annealing, at a lower temperature.

;;; The CHANGE procedure takes a possible solution and returns two values:
;;; (1) a procedure that, if applied to that solution, would modify it
;;; slightly, and (2) an exact real number indicating the difference
;;; between the value of the modified solution and the value of the
;;; original one.

;;; The COOL procedure takes an approximation and returns another such
;;; approximation.  The result should be less than the argument.  This
;;; procedure will be applied to decrease the simulated temperature (which
;;; is initially set to be much greater than the expected value of a
;;; proposed modification to the initial solution) from one annealing phase
;;; to the next.

(define (simulated-annealing
          initial-solution trial-limit change-limit change cool)

  ;; The initial temperature is computed by finding the greatest of the
  ;; values of TRIAL-LIMIT modifications that might be made to the initial
  ;; solution and doubling it.  This is supposed to ensure that even
  ;; modifications that would result in a large negative change in value
  ;; are occasionally approved during the initial annealing phases.

  (let ((initial-temperature
         (double (fold-natural trial-limit
                               (lambda (so-far)
                                 (receive (ignored delta)
                                          (change initial-solution)
                                   (max delta so-far)))
                               0)))

        ;; One step in the annealing process is completed when either the
        ;; number of proposed modifications or the number of modifications
        ;; actually performed reaches its upper bound

        (step-finished?
         (lambda (trials changes)
           (or (= trials trial-limit)
               (= changes change-limit))))

        ;; During one step in the annealing process, the temperature is
        ;; constant.  At that temperature, one cycle consists of generating
        ;; and evaluating a possible modification, submitting it to the
        ;; CHANGE-APPROVED? predicate to determine whether it should
        ;; actually be made, and making it if it is approved.

        (consider-alternative
         (lambda (temperature)
           (lambda (solution trials changes)
             (receive (changer delta) (change solution)
               (if (change-approved? delta temperature)
                   (values (changer solution)
                           (successor trials)
                           (successor changes)) 
                   (values solution (successor trials) changes)))))))

    ;; The ANNEALING-STEP procedure carries the process through all of the
    ;; proposed modifications that will be considered at one particular
    ;; temperature, then decreases the temperature for the next step.

    (let ((annealing-step
           (lambda (current temperature ignored)
             (let ((alternative (consider-alternative temperature)))
               (letrec ((trial
                         (lambda (solution trials changes)
                           (if (step-finished? trials changes)
                               (values solution
                                       (cool temperature)
                                       (zero? changes))
                               ((pipe alternative trial)
                                solution trials changes)))))
                 (trial current 0 0))))))

      ;; This body repeatedly invokes ANNEALING-STEP until some temperature
      ;; is reached at which none of the proposed modifications is
      ;; approved.

      ((until (project 2) annealing-step)
       initial-solution
       initial-temperature
       #f))))

;;; The CHANGE-APPROVED? procedure determines whether one possible solution
;;; should be replaced by another, given the amount by which the value of
;;; the new solution exceeds the value of the old one (DELTA) and the
;;; current simulated temperature.

;;; If the new solution's value is greater than or equal to the old one's
;;; (so that DELTA is non-negative), the modification is always approved.
;;; Otherwise, the probability that a modification is approved is the value
;;; of an exponential function with an exponent whose magnitude is directly
;;; proportional to DELTA (so that solutions that would greatly decrease
;;; the value are seldom approved) and inversely proportional to the
;;; simulated TEMPERATURE (so that new solutions are often approved at high
;;; temperatures, even when they decrease the value, but seldom approved at
;;; low temperatures).

(define (change-approved? delta temperature)
  (or (<= 0.0 delta)
      (< (random 1.0) (exp (/ delta temperature)))))

;;; In this example, we'll use simulated annealing to solve an instance of
;;; the traveling-salesperson problem.  The cities to be visited are the
;;; capitals of the nations of Africa; we'll represent each one as a tuple
;;; containing the city's name, latitude, and longitude.

(define (make-latitude degrees minutes direction)
  (list degrees minutes direction))

(define latitude-degrees (next-section list-ref 0))
(define latitude-minutes (next-section list-ref 1))
(define latitude-direction (next-section list-ref 2))

(define (latitude? something)
  (and (list? something)
       (= (length something) 3)
       (let ((first (list-ref something 0))
             (second (list-ref something 1))
             (third (list-ref something 2)))
         (and (natural-number? first)
              (<= first 90)
              (natural-number? second)
              (< second 60)
              (or (eq? third 'N) (eq? third 'S))))))

(define (make-longitude degrees minutes direction)
  (list degrees minutes direction))

(define longitude-degrees (next-section list-ref 0))
(define longitude-minutes (next-section list-ref 1))
(define longitude-direction (next-section list-ref 2))

(define (longitude? something)
  (and (list? something)
       (= (length something) 3)
       (let ((first (list-ref something 0))
             (second (list-ref something 1))
             (third (list-ref something 2)))
         (and (natural-number? first)
              (<= first 180)
              (natural-number? second)
              (< second 60)
              (or (eq? third 'E) (eq? third 'W))))))

(define (make-location name latitude longitude)
  (list name latitude longitude))

(define location-name (next-section list-ref 0))
(define location-latitude (next-section list-ref 1))
(define location-longitude (next-section list-ref 2))

(define (location? something)
  (and (list? something)
       (= (length something) 3)
       (and (string? (list-ref something 0))
            (latitude? (list-ref something 1))
            (longitude? (list-ref something 2)))))

;;; Here's our roster of African capitals.  (Two different nations named
;;; `Congo' are listed; the one with the capital in Kinshasa is the
;;; Democratic Republic of Congo.)

(define roster
  (vector
   (make-location "Abuja, Nigeria"
                  (make-latitude 9 12 'N)
                  (make-longitude 7 11 'E))
   (make-location "Accra, Ghana"
                  (make-latitude 5 33 'N)
                  (make-longitude 0 13 'W))
   (make-location "Addis Ababa, Ethiopia"
                  (make-latitude 9 0 'N)
                  (make-longitude 38 50 'E))
   (make-location "Algiers, Algeria"
                  (make-latitude 36 47 'N)
                  (make-longitude 3 3 'E))
   (make-location "Antananarivo, Madagascar"
                  (make-latitude 18 55 'S)
                  (make-longitude 47 31 'E))
   (make-location "Asmara, Eritrea"
                  (make-latitude 15 20 'N)
                  (make-longitude 38 53 'E))
   (make-location "Bamako, Mali"
                  (make-latitude 12 39 'N)
                  (make-longitude 8 0 'W))
   (make-location "Bangui, Central African Republic"
                  (make-latitude 4 22 'N)
                  (make-longitude 18 35 'E))
   (make-location "Banjul, Gambia"
                  (make-latitude 13 28 'N)
                  (make-longitude 16 39 'W))
   (make-location "Bissau, Guinea-Bissau"
                  (make-latitude 11 51 'N)
                  (make-longitude 15 35 'W))
   (make-location "Brazzaville, Congo"
                  (make-latitude 4 16 'S)
                  (make-longitude 15 17 'E))
   (make-location "Bujumbura, Burundi"
                  (make-latitude 3 22 'S)
                  (make-longitude 31 22 'E))
   (make-location "Cairo, Egypt"
                  (make-latitude 30 3 'N)
                  (make-longitude 31 15 'E))
   (make-location "Conakry, Guinea"
                  (make-latitude 9 31 'N)
                  (make-longitude 13 43 'W))
   (make-location "Dakar, Senegal"
                  (make-latitude 14 40 'N)
                  (make-longitude 17 26 'W))
   (make-location "Dar-es-Salaam, Tanzania"
                  (make-latitude 6 48 'S)
                  (make-longitude 39 17 'E))
   (make-location "Djibouti, Djibouti"
                  (make-latitude 11 36 'N)
                  (make-longitude 43 9 'E))
   (make-location "Freetown, Sierra Leone"
                  (make-latitude 8 30 'N)
                  (make-longitude 13 15 'W))
   (make-location "Gaborone, Botswana"
                  (make-latitude 24 45 'S)
                  (make-longitude 25 55 'E))
   (make-location "Harare, Zimbabwe"
                  (make-latitude 17 50 'S)
                  (make-longitude 31 3 'E))
   (make-location "Kampala, Uganda"
                  (make-latitude 0 19 'N)
                  (make-longitude 32 25 'E))
   (make-location "Khartoum, Sudan"
                  (make-latitude 15 36 'N)
                  (make-longitude 32 32 'E))
   (make-location "Kigali, Rwanda"
                  (make-latitude 1 57 'S)
                  (make-longitude 30 4 'E))
   (make-location "Kinshasa, Congo"
                  (make-latitude 4 18 'S)
                  (make-longitude 15 18 'E))
   (make-location "Libreville, Gabon"
                  (make-latitude 0 23 'N)
                  (make-longitude 9 27 'E))
   (make-location "Lilongwe, Malawi"
                  (make-latitude 13 59 'S)
                  (make-longitude 33 44 'E))
   (make-location "Lome, Togo"
                  (make-latitude 6 8 'N)
                  (make-longitude 1 13 'E))
   (make-location "Luanda, Angola"
                  (make-latitude 8 48 'S)
                  (make-longitude 13 14 'E))
   (make-location "Lusaka, Zambia"
                  (make-latitude 15 25 'S)
                  (make-longitude 28 17 'E))
   (make-location "Malabo, Equatorial Guinea"
                  (make-latitude 3 45 'N)
                  (make-longitude 8 47 'E))
   (make-location "Maputo, Mozambique"
                  (make-latitude 25 58 'S)
                  (make-longitude 32 35 'E))
   (make-location "Maseru, Lesotho"
                  (make-latitude 29 28 'S)
                  (make-longitude 27 30 'E))
   (make-location "Mbabane, Swaziland"
                  (make-latitude 26 18 'S)
                  (make-longitude 31 6 'E))
   (make-location "Mogadishu, Somalia"
                  (make-latitude 2 1 'N)
                  (make-longitude 45 20 'E))
   (make-location "Monrovia, Liberia"
                  (make-latitude 6 18 'N)
                  (make-longitude 10 47 'W))
   (make-location "Moroni, Comoros"
                  (make-latitude 11 41 'S)
                  (make-longitude 43 16 'E))
   (make-location "Nairobi, Kenya"
                  (make-latitude 1 17 'S)
                  (make-longitude 36 49 'E))
   (make-location "Ndjamena, Chad"
                  (make-latitude 12 7 'N)
                  (make-longitude 15 3 'E))
   (make-location "Niamey, Niger"
                  (make-latitude 13 31 'N)
                  (make-longitude 2 7 'E))
   (make-location "Nouakchott, Mauritania"
                  (make-latitude 18 6 'N)
                  (make-longitude 15 57 'W))
   (make-location "Ouagadougou, Burkina Faso"
                  (make-latitude 12 22 'N)
                  (make-longitude 1 31 'W))
   (make-location "Port Louis, Mauritius"
                  (make-latitude 20 10 'S)
                  (make-longitude 57 30 'E))
   (make-location "Porto-Novo, Benin"
                  (make-latitude 6 29 'N)
                  (make-longitude 2 37 'E))
   (make-location "Praia, Cape Verde"
                  (make-latitude 14 55 'N)
                  (make-longitude 23 31 'W))
   (make-location "Pretoria, South Africa"
                  (make-latitude 25 45 'S)
                  (make-longitude 28 10 'E))
   (make-location "Rabat, Morocco"
                  (make-latitude 34 2 'N)
                  (make-longitude 6 51 'W))
   (make-location "Sao Tome, Sao Tome and Principe"
                  (make-latitude 0 20 'N)
                  (make-longitude 6 44 'E))
   (make-location "Tripoli, Libya"
                  (make-latitude 32 54 'N)
                  (make-longitude 15 11 'E))
   (make-location "Tunis, Tunisia"
                  (make-latitude 36 48 'N)
                  (make-longitude 10 11 'E))
   (make-location "Victoria, Seychelles"
                  (make-latitude 4 38 'S)
                  (make-longitude 55 27 'E))
   (make-location "Windhoek, Namibia"
                  (make-latitude 22 34 'S)
                  (make-longitude 17 6 'E))
   (make-location "Yamoussoukro, Cote d'Ivoire"
                  (make-latitude 6 49 'N)
                  (make-longitude 5 17 'W))
   (make-location "Yaounde, Cameroon"
                  (make-latitude 3 52 'N)
                  (make-longitude 11 31 'E))))

;;; We'll often need the number of cities in this list.

(define roster-length (vector-length roster))

;;; In this example, we'll assume that the cost of travelling from one city
;;; to another is proportional to the great-circle distance between them,
;;; as computed from their latitudes and longitudes.  The
;;; GREAT-CIRCLE-DISTANCE procedure carries out that computation, returning
;;; distances in kilometers.  Since trigonometric functions return
;;; approximations rather than numbers, the initially computed distance is
;;; also an approximation.  GREAT-CIRCLE-DISTANCE rounds it to the nearest
;;; kilometer and restores exactness to the result.

(define earth-radius 6367)  ; the earth's radius, in kilometers

(define (great-circle-distance unua dua)
  (let ((unua-latitude (latitude->radians (location-latitude unua)))
        (unua-longitude (longitude->radians (location-longitude unua)))
        (dua-latitude (latitude->radians (location-latitude dua)))
        (dua-longitude (longitude->radians (location-longitude dua))))
    (let ((a (+ (square (sin (/ (- dua-latitude unua-latitude) 2)))
                (* (cos unua-latitude)
                   (cos dua-latitude)
                   (square (sin (/ (- dua-longitude unua-longitude)
                                   2)))))))
      (inexact->exact
        (round (* earth-radius 2 (atan (sqrt a) (sqrt (- 1 a)))))))))

;;; Latitudes are given in degrees and minutes north or south of the
;;; equator, and longitude in degrees and minutes east or west of the
;;; Greenwich meridian.  The following procedures convert them to radians.

(define pi 3.141592653589793)

(define (latitude->radians lat)
  (let ((magnitude
         (* (+ (latitude-degrees lat) (/ (latitude-minutes lat) 60))
            (/ pi 180))))
    (if (eq? (latitude-direction lat) 'N)
        magnitude
        (- magnitude))))

(define (longitude->radians long)
  (let ((magnitude
         (* (+ (longitude-degrees long) (/ (longitude-minutes long) 60))
            (/ pi 180))))
    (if (eq? (longitude-direction long) 'E)
        magnitude
        (- magnitude))))

;;; Instead of recomputing the distances repeatedly, we'll compute each one
;;; once and store the result in a two-dimensional array.

(define distance-array
  (generate-vector
    (lambda (origin-index)
      (let ((origin (vector-ref roster origin-index)))
        (generate-vector
          (lambda (destination-index)
            (let ((destination (vector-ref roster destination-index)))
              (great-circle-distance origin destination)))
          roster-length)))
    roster-length))

;;; Because we'll be accessing this particular array repeatedly, it's
;;; convenient to define a special procedure DISTANCE that takes the serial
;;; numbers of two cities (i.e., their positions in the roster) and returns
;;; the distance between them, in kilometers.

(define (distance origin-index destination-index)
  (vector-ref (vector-ref distance-array origin-index) destination-index))

;;; Now we're almost ready to implement the CHANGE, COOL, and
;;; INITIAL-SOLUTION procedures that characterize our problem.  We'll
;;; represent each solution as a pair in which the car is a possible tour
;;; of the cities on the roster and the cdr is the value of that tour.

(define make-solution cons)
(define solution-tour car)
(define solution-value cdr)

;;; A tour can be represented as a vector, equal in length to ROSTER, in
;;; which each entry is the position in ROSTER of a different city; we'll
;;; assume that the first trip is from the city identified by the first
;;; element of the tour to the city identified by the second, the second
;;; trip from the city identified by the second element of the tour to the
;;; city identified by the third, and so on, finishing with a trip from the
;;; city identified by the last element of the tour back to the city
;;; identified by the first.

;;; Because we'll often be treating this vector as a circular structure, it
;;; is helpful to define a procedure that maps any integer into the vector
;;; index that would result from going ``around the corner'' from one end
;;; of the vector to the other as many times as necessary to bring the
;;; index into range:

(define cyclic-index (next-section modulo roster-length))

;;; To construct the initial solution, we'll just imagine travelling to the
;;; cities in the order in which they appear in the roster above.  The
;;; value of this solution may be defined as the negative of the sum of the
;;; costs of the individual city-to-city trips that make it up.  (We're
;;; trying to make this negated sum as great -- as nearly positive -- as
;;; possible, by making its magnitude as small as possible.)

(define initial-solution
  (let ((tour (generate-vector identity roster-length)))
    (make-solution tour
                   (- (lower-ply-natural
                        (predecessor roster-length)
                        (lambda (index sum)
                          (+ sum
                             (distance (vector-ref tour index)
                                       (vector-ref tour (successor index)))))
                        (distance (vector-ref tour
                                              (predecessor roster-length))
                                  (vector-ref tour 0)))))))

;;; We'll use the standard recommended cooling procedure.

(define cool (initial-section * 0.9))

;;; We'll consider two kinds of modifications: (1) reversal, in which we
;;; break out a segment of the tour, reverse it end for end, and then
;;; splice it back in at the original position, and (2) transposal, in
;;; which we again break out a segment of the tour, splice over the break,
;;; open up the tour at some other point between two cities, and insert the
;;; segment that we previously removed.

;;; If we start with a tour like A-B-C-D-E-F-G-H-I-A, then breaking at C
;;; and G and reversing the resulting segment C-D-E-F would yield
;;; A-B-F-E-D-C-G-H-I-A.  Breaking at B and E and transposing the segment
;;; B-C-D to a break at H would yield A-E-F-G-B-C-D-H-I-A.

;;; CHANGE tosses a coin to decide whether to propose a reversal or a
;;; transposal.

(define (change solution)
  ((if (zero? (random 2)) reversal transposal) solution))

;;; REVERSAL selects one endpoint for the segment to be reversed and
;;; chooses the length of the segment, being careful to create a segment
;;; containing at least two cities and to leave at least two cities in the
;;; unreversed portion of the tour.  It then finds the other endpoint.

;;; Once the segment to be reversed is identified, it does not matter
;;; whether we reverse it or reverse the rest of the tour, since the new
;;; tours that result will differ only the direction in which each trip is
;;; taken.  Hence we lose no generality in requiring that START, the
;;; beginning (inclusive) endpoint of the segment, is less than FINISH, the
;;; ending (exclusive) endpoint.

;;; The amount by which the value of the new tour exceeds the value of the
;;; old is computed by adding the distances of the city-to-city trips
;;; removed from the tour to the negatives of the distances of the trips
;;; added to the tour.  (If the signs seem backwards, remember that the
;;; value of a tour is the _negative_ of the sum of the distances of the
;;; trips that it comprises.)

(define (reversal solution)
  (let ((tour (solution-tour solution))
        (one-end (random roster-length))
        (reversed-segment-length
         (+ (random (- roster-length 3)) 2)))
    (let ((other-end (cyclic-index (+ one-end reversed-segment-length))))
      (let ((start (min one-end other-end))
            (finish (max one-end other-end)))
        (let ((pre-start (cyclic-index (predecessor start)))
              (pre-finish (cyclic-index (predecessor finish))))
          (let ((delta (+ (distance (vector-ref tour pre-start)
                                    (vector-ref tour start))
                          (distance (vector-ref tour pre-finish)
                                    (vector-ref tour finish))
                          (- (distance (vector-ref tour pre-start)
                                       (vector-ref tour pre-finish)))
                          (- (distance (vector-ref tour start)
                                       (vector-ref tour finish))))))
            (values (reverser start finish delta) delta)))))))

;;; REVERSER takes as arguments the indices of the (inclusive) starting and
;;; (exclusive) ending points of a tour segment to be reversed, and an
;;; indication of how making the reversal affects the value of the tour.
;;; It returns a procedure that, when applied to a solution, will construct
;;; a new solution, with the specified segment reversed and the value
;;; adjusted as indicated.

(define (reverser start finish delta)
  (lambda (solution)
    (let ((old-tour (solution-tour solution))
          (old-value (solution-value solution)))
      (make-solution
        (generate-vector (lambda (index)
                           (vector-ref old-tour
                                       (if (and (<= start index)
                                                (< index finish))
                                           (+ start (- finish index 1))
                                           index)))
                         roster-length)
        (+ old-value delta)))))

;;; TRANSPOSAL selects a starting point for the segment to be transposed
;;; and chooses the length of the segment, being careful to create a
;;; segment containing at least one city and to leave at least two cities
;;; in the untransposed portion of the tour.  It then finds the (exclusive)
;;; finish point of the segment and randomly selects a ``break'' point, not
;;; within the segment, to open up the rest of the tour for the insertion
;;; of the transposed segment.

;;; The amount by which the value of the new tour exceeds the value of the
;;; old is computed, as above, by adding the distances of the city-to-city
;;; trips removed from the tour to the negatives of the distances of the
;;; trips added to the tour.  In this case, there are three of each -- the
;;; trips ending at the cities in positions START, FINISH, and BREAK will
;;; all begin in different places after the transposition has been made.

(define (transposal solution)
  (let ((tour (solution-tour solution))
        (start (random roster-length))
        (splice-length (successor (random (- roster-length 2)))))
    (let ((finish (cyclic-index (+ start splice-length)))
          (break
           (cyclic-index (- start (successor
                                    (random (predecessor
                                              (- roster-length
                                                 splice-length))))))))
      (let ((pre-start (cyclic-index (predecessor start)))
            (pre-finish (cyclic-index (predecessor finish)))
            (pre-break (cyclic-index (predecessor break))))
        (let ((delta (+ (distance (vector-ref tour pre-start)
                                  (vector-ref tour start))
                        (distance (vector-ref tour pre-finish)
                                  (vector-ref tour finish))
                        (distance (vector-ref tour pre-break)
                                  (vector-ref tour break))
                        (- (distance (vector-ref tour pre-start)
                                     (vector-ref tour finish)))
                        (- (distance (vector-ref tour pre-finish)
                                     (vector-ref tour break)))
                        (- (distance (vector-ref tour pre-break)
                                     (vector-ref tour start))))))
          (values (transposer start finish break delta) delta))))))

;;; TRANSPOSER takes as arguments the indices of the (inclusive) starting
;;; and (exclusive) ending points of a tour segment to be transposed, the
;;; index of the position at which the segment is to be reinserted, and an
;;; indication of how making the transposal affects the value of the tour.
;;; It returns a procedure that, when applied to a solution, will construct
;;; a new solution, with the specified segment transposed to the new
;;; location and the value adjusted as indicated.

;;; At first it might seem necessary to write different code for the three
;;; possible orderings of START, FINISH, and BREAK:

;;; Case 1: Change
;;;
;;; 0 ---| START --------| FINISH -----| BREAK ------------| ROSTER-LENGTH
;;;
;;; to
;;;
;;; 0 ---| FINISH -----| START --------| BREAK ------------| ROSTER-LENGTH

;;; Case 2: Change
;;;
;;; 0 ---| BREAK ------------| START --------| FINISH -----| ROSTER-LENGTH
;;;
;;; to
;;;
;;; 0 ---| START --------| BREAK ------------| FINISH -----| ROSTER-LENGTH

;;; Case 3 (``around the corner''): Change
;;;
;;; 0 ---| FINISH -----| BREAK ------------| START --------| ROSTER-LENGTH
;;;
;;; to
;;;
;;; 0 ---| BREAK ------------| FINISH -----| START --------| ROSTER-LENGTH

;;; However, the symmetry of these cases makes it possible to factor out a
;;; common transposal method (here called ROTOR), which we can provide with
;;; the appropriate parameters to obtain the particular transposer we want.

(define (transposer start finish break delta)
  (lambda (solution)
    (let ((old-tour (solution-tour solution))
          (old-value (solution-value solution)))
      (let ((rotor
             (lambda (fore mid aft)
               (lambda (index)
                 (vector-ref old-tour
                             (if (< index fore)
                                 index
                                 (if (< index (+ fore (- aft mid)))
                                     (+ mid (- index fore))
                                     (if (< index aft)
                                         (+ fore (- index
                                                    (+ fore (- aft mid))))
                                         index))))))))
        (make-solution (generate-vector
                         (if (< start finish break)
                             (rotor start finish break)
                             (if (< break start finish)
                                 (rotor break start finish)
                                 (if (< finish break start)
                                     (rotor finish break start)
                                     (error "This can't happen!"))))
                         roster-length)
                       (+ old-value delta))))))

;;; Once we have a solution, we'll want to print it out.  The DISPLAY-TRIP
;;; procedure, which we'll invoke once for each trip in the tour, prints
;;; out the names of the cities of origin and destination, given their
;;; indices in ROSTER, and the distance between them in kilometers, as
;;; stored in DISTANCE-ARRAY.

(define (display-trip origin-index destination-index)
  (let ((origin (vector-ref roster origin-index))
        (destination (vector-ref roster destination-index)))
    (display "From ")
    (display (location-name origin))
    (display " to ")
    (display (location-name destination))
    (display " (")
    (display (distance origin-index destination-index))
    (display " km)")
    (newline)))

;;; To display the final solution, we print out each trip in the tour and
;;; finally give the total distance travelled.

(define (display-solution solution)
  (let ((tour (solution-tour solution))
        (value (solution-value solution)))
    (lower-ply-natural (predecessor roster-length)
                       (lambda (index)
                         (display-trip (vector-ref tour index)
                                       (vector-ref tour (successor index)))
                         (black-hole)))
    (display-trip (vector-ref tour (predecessor roster-length))
                  (vector-ref tour 0))
    (newline)
    (display "length of tour = ")
    (display (- value))
    (display " km")
    (newline)))

;;; We now invoke SIMULATED-ANNEALING with parameters, chosen somewhat
;;; arbitrarily.  TRIAL-LIMIT is set to five hundred times the number of
;;; cities to be visited, CHANGE-LIMIT to ten times that number.

(display-solution (simulated-annealing initial-solution
                                       (* 500 roster-length)
                                       (* 10 roster-length)
                                       change
                                       cool))

;;; A test run produced the following output:

;;; From Malabo, Equatorial Guinea to Sao Tome, Sao Tome and Principe (443 km)
;;; From Sao Tome, Sao Tome and Principe to Libreville, Gabon (302 km)
;;; From Libreville, Gabon to Brazzaville, Congo (829 km)
;;; From Brazzaville, Congo to Kinshasa, Congo (4 km)
;;; From Kinshasa, Congo to Luanda, Angola (550 km)
;;; From Luanda, Angola to Windhoek, Namibia (1584 km)
;;; From Windhoek, Namibia to Gaborone, Botswana (929 km)
;;; From Gaborone, Botswana to Pretoria, South Africa (252 km)
;;; From Pretoria, South Africa to Maseru, Lesotho (418 km)
;;; From Maseru, Lesotho to Mbabane, Swaziland (499 km)
;;; From Mbabane, Swaziland to Maputo, Mozambique (153 km)
;;; From Maputo, Mozambique to Harare, Zimbabwe (918 km)
;;; From Harare, Zimbabwe to Lusaka, Zambia (399 km)
;;; From Lusaka, Zambia to Lilongwe, Malawi (607 km)
;;; From Lilongwe, Malawi to Bujumbura, Burundi (1208 km)
;;; From Bujumbura, Burundi to Kigali, Rwanda (214 km)
;;; From Kigali, Rwanda to Kampala, Uganda (363 km)
;;; From Kampala, Uganda to Nairobi, Kenya (520 km)
;;; From Nairobi, Kenya to Dar-es-Salaam, Tanzania (671 km)
;;; From Dar-es-Salaam, Tanzania to Moroni, Comoros (697 km)
;;; From Moroni, Comoros to Antananarivo, Madagascar (924 km)
;;; From Antananarivo, Madagascar to Port Louis, Mauritius (1055 km)
;;; From Port Louis, Mauritius to Victoria, Seychelles (1740 km)
;;; From Victoria, Seychelles to Mogadishu, Somalia (1345 km)
;;; From Mogadishu, Somalia to Addis Ababa, Ethiopia (1058 km)
;;; From Addis Ababa, Ethiopia to Djibouti, Djibouti (553 km)
;;; From Djibouti, Djibouti to Asmara, Eritrea (620 km)
;;; From Asmara, Eritrea to Khartoum, Sudan (681 km)
;;; From Khartoum, Sudan to Cairo, Egypt (1611 km)
;;; From Cairo, Egypt to Tripoli, Libya (1554 km)
;;; From Tripoli, Libya to Tunis, Tunisia (629 km)
;;; From Tunis, Tunisia to Algiers, Algeria (635 km)
;;; From Algiers, Algeria to Rabat, Morocco (947 km)
;;; From Rabat, Morocco to Nouakchott, Mauritania (1988 km)
;;; From Nouakchott, Mauritania to Praia, Cape Verde (880 km)
;;; From Praia, Cape Verde to Dakar, Senegal (654 km)
;;; From Dakar, Senegal to Banjul, Gambia (158 km)
;;; From Banjul, Gambia to Bissau, Guinea-Bissau (214 km)
;;; From Bissau, Guinea-Bissau to Conakry, Guinea (330 km)
;;; From Conakry, Guinea to Freetown, Sierra Leone (124 km)
;;; From Freetown, Sierra Leone to Monrovia, Liberia (366 km)
;;; From Monrovia, Liberia to Yamoussoukro, Cote d'Ivoire (610 km)
;;; From Yamoussoukro, Cote d'Ivoire to Bamako, Mali (713 km)
;;; From Bamako, Mali to Ouagadougou, Burkina Faso (704 km)
;;; From Ouagadougou, Burkina Faso to Niamey, Niger (414 km)
;;; From Niamey, Niger to Accra, Ghana (921 km)
;;; From Accra, Ghana to Lome, Togo (171 km)
;;; From Lome, Togo to Porto-Novo, Benin (159 km)
;;; From Porto-Novo, Benin to Abuja, Nigeria (586 km)
;;; From Abuja, Nigeria to Ndjamena, Chad (918 km)
;;; From Ndjamena, Chad to Bangui, Central African Republic (945 km)
;;; From Bangui, Central African Republic to Yaounde, Cameroon (785 km)
;;; From Yaounde, Cameroon to Malabo, Equatorial Guinea (303 km)
;;; 
;;; length of tour = 36855 km