{ traveler

  A traveling salesman is planning a tour of n cities in which he has
  business.  He wishes to visit each city exactly once, to go directly from
  one city to another, and to end the tour in the same city that he starts
  out from.  He has a table showing the cost of traveling from any city on
  his itinerary to any other.  He wants to make the tour for the smallest
  possible total cost.  How should he determine the best route?

  One possibility would be to generate each of the n! possible routes,
  computing the cost of each one, and to select the minimum.  However, if n
  is not very small, this strategy is computationally infeasible.

  The method used in this program is to construct an initial path at
  random, and then to perturb it in various random ways, accepting all
  perturbations that decrease the total cost and rejecting some that increase
  it, particularly if the increase is large.  (The reason for accepting some
  cost-increasing changes is to make it less likely that the algorithm will
  get trapped in a local minimum that is not close to a global minimum.)  The
  algorithm calls for a progressively stronger barrier against
  cost-increasing perturbations, until it becomes unlikely that any further
  perturbations will lead to improvements in the path.

  A fuller description of the algorithm appears in section 10.9 of
  _Numerical recipes_, by William H. Press et al. (Cambridge: Cambridge
  University Press, 1986), pages 326-34, as the "Method of Simulated
  Annealing."

  Programmer: John Stone, Grinnell College
  Modification history:
      Original version -- April 28, 1986
      Revised for Suns -- May 9, 1988
      More accurate distance function introduced -- April 26, 1990
}
program traveler (input, output);
const
    MAXVERTICES = 200;
        { the largest number of cities that can figure in the salesman's 
          itinerary }
    MAXSTRING = 30;
        { maximum number of characters in the name of a city }
type
    VERTEXCOUNT = 0 .. MAXVERTICES;
        { the number of vertices in a graph }
    STRRANGE = 1 .. MAXSTRING;
        { positions within a string }
    STR = packed array [STRRANGE] of char;
        { a standard Pascal string }
    STRCOUNT = 0 .. MAXSTRING;
        { possible lengths of a string }
    NAME = record
               chars: STR;
               len: STRCOUNT
           end;
        { the name of a city, with an indication of the length of the name }
    VERTEXRANGE = 1 .. MAXVERTICES;
        { used for enumerating vertices in a graph }
    VERTEX = record
                 vname: NAME;
                 latitude: real;
                 longitude: real
             end;
        { the coordinates and identification of a city on the tour }
    GRAPH = record
                items: array [VERTEXRANGE] of VERTEX;
                size: VERTEXCOUNT
            end;
        { a collection of vertices, within which an itinerary is to be
          found }
    COSTTABLE = array [VERTEXRANGE, VERTEXRANGE] of real;
        { a table giving costs (= distances = edge weights) of traveling
          between cities (= vertices) }
    PATHLIST = array [VERTEXRANGE] of VERTEXRANGE;
        { an object of this type is a list of the serial numbers of the
          vertices along a tentative itinerary, subject to revision, in
          the order in which they are to be visited }
    PATH = record
               vertices: VERTEXCOUNT;
               tour: PATHLIST;
               cost: real
           end;
        { an itinerary and its total cost }
var
    map: GRAPH;
        { information about the cities to be covered }
    distance: COSTTABLE;
        { the costs of traveling between cities }
    current: PATH;
        { a tentative solution to be developed further }
    maxtrials: integer;
        { the maximum number of attempts to change the path that will be 
          made in each phase of the simulation }
    maxchanges: integer;
        { the maximum number of changes in the path that will be accepted 
          in each phase of the simulation }
    mutability: real;
        { a measure of the extent to which path changes that increase the 
          total cost will be permitted }
    trials: integer;
        { counts off the number of attempts to modify the path at a given 
          level of mutability }
    changemade: Boolean;
        { indicates whether a particular trial resulted in a modification 
          to the path }
    changes: integer;
        { the number of changes made at a given level of mutability }

    { readstring

      This procedure collects one string from a line of standard input and 
      measures its length.
    }
    procedure readstring (var s: NAME);
    var
        ch: char;
            { one character at a time from standard input }
    begin
        with s do begin
            len := 0;
            while not eoln and (len < MAXSTRING) do begin
                len := len + 1;
                read (ch);
                chars[len] := ch
            end;
            readln
        end
    end;

    { readcoord

      This procedure reads in a latitude or longitude (in the form
      
                            degrees,minutes

      where degrees is an integer between -180 and 180 and minutes is an
      integer between 0 and 59), converts it to radians, and returns the
      result through its parameter.
    }
    procedure readcoord (var coord: real);
    const
        PI = 3.14159265358979324;
    var
        degrees, minutes: integer;
        separator: char;
            { the comma between the degrees and the minutes }
    begin
        readln (degrees, separator, minutes);
        if degrees < 0 then
            minutes := -minutes;
        coord := (PI / 180.0) * (degrees + minutes / 60.0)
    end;

    { readproblem

      This procedure collects all the information needed to set up the
      problem.
    }
    procedure readproblem (var problem: GRAPH);
    begin
        with problem do begin
            size := 0;
            while not eof and (size < MAXVERTICES) do begin
                size := size + 1;
                with items[size] do begin
                    readstring (vname);
                    readcoord (latitude);
                    readcoord (longitude)
                end
            end
        end
    end;

    { arccos

      This function computes and returns the arc cosine, in radians, of its
      argument, which must be a real number between -1 and 1 inclusive.
    }
    function arccos (x: real): real;
    const
        EPSILON = 1.0e-20;
            { a real number that is effectively equal to zero; if x differs
              from 1 or from -1 by less than EPSILON, no attempt is made to
              compute the value of the arc cosine by the standard formula,
              for fear of encountering a division by zero }
        PI = 3.14159265358979324;
    begin
        if abs (x - 1.0) < EPSILON then
            arccos := 0.0
        else if abs (x + 1.0) < EPSILON then
            arccos := PI
        else
            arccos := (PI / 2.0) - arctan (x / sqrt (1.0 - sqr (x)))
    end;

    { computedistances

      This procedure fills in the table that gives the cost of traveling
      from any vertex to any other.
    }
    procedure computedistances (var map: GRAPH; var distance: COSTTABLE);
    const
        RADIUS = 6.377e3;
            { the radius of the earth, in kilometers }
    var
        start, finish: VERTEXRANGE;
            { two vertices whose separation is to be computed }
    begin
        with map do
            for start := 1 to size do begin
                distance[start, start] := 0.0;
                for finish := start + 1 to size do begin
                    distance[start, finish] :=
                        RADIUS * arccos (sin (items[start].latitude) *
                                         sin (items[finish].latitude) +
                                         cos (items[start].latitude) *
                                         cos (items[finish].latitude) *
                                         cos (items[finish].longitude -
                                                  items[start].longitude));
                    distance[finish, start] := distance[start, finish]
                end
            end
    end;

    { initialpath

      This procedure sets up an initial path in the most straightforward
      way possible.
    }
    procedure initialpath (var map: GRAPH; var distance: COSTTABLE;
        var p: PATH);
    var
        city: VERTEXRANGE;
            { counts off the cities to be visited }
    begin
        with map, p do begin
            vertices := size;
            cost := distance[size, 1];
            for city := 1 to size do begin
                tour[city] := city;
                cost := cost + distance[city, city + 1]
            end
        end
    end;
            
    { randomint

      This function returns a random integer in a specified range.
    }
    function randomint (lower, upper: integer): integer;
    begin
        randomint := lower + trunc (random(0) * (upper - lower + 1))
    end;

    { randomize

      This procedure randomly sets the seed for the random-number generator.
    }
    procedure randomize;
    var
        dummy: integer;
    begin
        dummy := seed (wallclock)
    end;

    { reversal

      This procedure randomly chooses a segment of the given path and
      determines the change that would occur in the total cost of the path
      if the segment were spliced out, reversed, and spliced back in again.
      This change is returned as the value of the function, and the two ends
      of the segment that would be reversed are returned through the
      parameters start and finish.
    }
    function reversal (var p: PATH; var distance: COSTTABLE;
        var start, finish: VERTEXRANGE): real;
    var
        length: VERTEXRANGE;
            { the number of edges in the path segment to be reversed -- 
              between 1 and vertices div 2 because longer reversals just 
              produce the reverse paths of shorter ones }
        prestart: VERTEXRANGE;
            { the number of the vertex just before 'start' along the path }
        postfinish: VERTEXRANGE;
            { the number of the vertex just after 'finish' along the path }
    begin
        with p do begin
            start := randomint (1, vertices);
            length := randomint (1, vertices div 2);
            if start + length <= vertices then
                finish := start + length
            else
                finish := start + length - vertices;
            if start = 1 then
                prestart := vertices
            else
                prestart := start - 1;
            if finish = vertices then
                postfinish := 1
            else
                postfinish := finish + 1;
            reversal := distance[tour[prestart], tour[finish]]
                      + distance[tour[start], tour[postfinish]]
                      - distance[tour[prestart], tour[start]]
                      - distance[tour[finish], tour[postfinish]]
        end
    end;

    { transposal

      This procedure randomly selects a segment of the given path and a
      vertex not contained in that segment and determines the change that
      would occur in the total cost of the path if the segment were spliced
      out from its current position and spliced in again between the other
      chosen vertex and its successor along the path.  This change is
      returned as the value of the function, and the two ends of the
      segment that would be transposed and the vertex after which it
      would be inserted are returned through the parameters start,
      finish, and break.
    }
    function transposal (var p: PATH; var distance: COSTTABLE;
        var start, finish, break: VERTEXRANGE): real;
    var
        length: VERTEXRANGE;
            { the length of the segment to be spliced in; must be no 
              greater than vertices - 3 in order to leave a path to splice 
              into }
        distancetobreak: VERTEXRANGE;
            { the number of positions between 'finish' and the randomly
              selected break point along the path }
        prestart: VERTEXRANGE;
            { the number of the vertex just before 'start' along the path }
        postfinish: VERTEXRANGE;
            { the number of the vertex just after 'finish' along the path }
        postbreak: VERTEXRANGE;
            { the number of the vertex just after 'break' along the path; 
              this is the vertex before which the splice will be made }
    begin
        with p do begin
            start := randomint (1, vertices);
            length := randomint (1, vertices - 3);
            if start + length <= vertices then
                finish := start + length
            else
                finish := start + length - vertices;
            distancetobreak := randomint (1, vertices - length - 2);
            if finish + distancetobreak <= vertices then
                break := finish + distancetobreak
            else
                break := finish + distancetobreak - vertices;
            if start = 1 then
                prestart := vertices
            else
                prestart := start - 1;
            if finish = vertices then
                postfinish := 1
            else
                postfinish := finish + 1;
            if break = vertices then
                postbreak := 1
            else
                postbreak := break + 1;
            transposal := distance[tour[prestart], tour[postfinish]]
                        + distance[tour[break], tour[start]]
                        + distance[tour[finish], tour[postbreak]]
                        - distance[tour[prestart], tour[start]]
                        - distance[tour[finish], tour[postfinish]]
                        - distance[tour[break], tour[postbreak]]
        end
    end;
        
    { findrange

      This function generates a number that should be safely larger than the 
      cost increment of any random reversal or transposal of a path segment.
    }
    function findrange (var p: PATH; var distance: COSTTABLE;
        trials: integer): real;
    var
        maxcost: real;
            { the largest cost increment encountered in 'trials' random 
              reversals and the same number of random transposals }
        trial: integer;
            { counts off the random reversals and transposals }
        cost: real;
            { the cost increment of one random reversal or transposal }
        start, finish, break: VERTEXRANGE;
            { unused values returned by reversal and transposal procedures }
    begin
        maxcost := 0.0;
        for trial := 1 to trials do begin
            cost := reversal (p, distance, start, finish);
            if cost > maxcost then
                maxcost := cost
        end;
        for trial := 1 to trials do begin
            cost := transposal (p, distance, start, finish, break);
            if cost > maxcost then
                maxcost := cost
        end;
        findrange := 2.0 * maxcost
    end;

    { Metropolis

      This function determines whether it is appropriate to perform a path
      mutation, given the size of the change that it would make in the
      total cost and the current level of permitted mutability.  Changes
      that decrease the total cost are always approved.  The Metropolis
      algorithm is used to decide whether a cost-increasing change is to
      be permitted:  It generates a random number in the range [0, 1) and
      compares it with the result of raising e to a fractional power, the
      numerator of the fraction being the negative of the change in cost (so
      that big cost increases are seldom approved) and the denominator being
      the mutability (so that very few changes are approved if the
      mutability is low).
    }
    function Metropolis (deltacost, mutability: real): Boolean;
    const
        BIGRATIO = 200;
            { A change that increases the cost by more than BIGRATIO
              times the current mutability should never be approved. }
    begin
        if deltacost <= 0.0 then
            Metropolis := TRUE
        else if BIGRATIO < deltacost / mutability then
            Metropolis := FALSE
        else
            Metropolis := (random(0) < exp (-deltacost / mutability))
    end;

    { reverse

      This procedure actually effects the reversal of a specified segment
      of a path.
    }
    procedure reverse (var p: PATH; start, finish: VERTEXRANGE;
                       deltacost: real);
    var
        v: integer;
            { counts off the points along the path }
        temp: PATHLIST;
            { temporary storage for the rearranged path, as it is 
              constructed }
        source: integer;
            { counts off points along the path in its original form as they 
              are transferred to different positions along the new path }
    begin
        with p do begin
            cost := cost + deltacost;
            if start < finish then begin
                for v := 1 to start - 1 do
                    temp[v] := tour[v];
                for v := start to finish do
                    temp[v] := tour[finish - v + start];
                for v := finish + 1 to vertices do
                    temp[v] := tour[v]
            end
            else begin
                source := start;
                for v := finish downto 1 do begin
                    temp[v] := tour[source];
                    if source = vertices then
                        source := 1
                    else
                        source := source + 1
                end;
                for v := vertices downto start do begin
                    temp[v] := tour[source];
                    if source = vertices then
                        source := 1
                    else
                        source := source + 1
                end;
                for v := finish + 1 to start - 1 do
                    temp[v] := tour[v]
            end;
            tour := temp
        end
    end;

    { splice

      This procedure actually effects a transposition operation in which
      the segment of path p lying between start and finish, inclusive,
      is removed from its present position and spliced back into the
      remainder of the path, starting immediately after break.
    }
    procedure splice (var p: PATH; start, finish, break: VERTEXRANGE;
                         deltacost: real);
    var
        v: integer;
            { counts off the points along the path }
        temp: PATHLIST;
            { temporary storage for the rearranged path, as it is 
              constructed }
        source: integer;
            { counts off points along the path in its original form as they 
              are transferred to different positions along the new path }
    begin
        with p do begin
            cost := cost + deltacost;
            if start < break then begin
                for v := 1 to start - 1 do
                    temp[v] := tour[v];
                v := start;
                for source := finish + 1 to break do begin
                    temp[v] := tour[source];
                    v := v + 1
                end;
                for source := start to finish do begin
                    temp[v] := tour[source];
                    v := v + 1
                end;
                for v := break + 1 to vertices do
                    temp[v] := tour[v]
            end
            else if break < finish then begin
                for v := 1 to break do
                    temp[v] := tour[v];
                v := break + 1;
                for source := start to finish do begin
                    temp[v] := tour[source];
                    v := v + 1
                end;
                for source := break + 1 to start - 1 do begin
                    temp[v] := tour[source];
                    v := v + 1
                end;
                for v := finish + 1 to vertices do
                    temp[v] := tour[v]
            end
            else begin
                for v := 1 to finish do
                    temp[v] := tour[v];
                v := finish + 1;
                for source := break + 1 to start - 1 do begin
                    temp[v] := tour[source];
                    v := v + 1
                end;
                for source := finish + 1 to break do begin
                    temp[v] := tour[source];
                    v := v + 1
                end;
                for v := start to vertices do
                    temp[v] := tour[v]
            end;
            tour := temp
        end
    end;

    { mutate

      This procedure selects a random reversal or transposal and, if the
      Metropolis algorithm approves, effects it.
    }
    procedure mutate (var p: PATH; mutability: real; var distance: COSTTABLE;
        var changemade: Boolean);
    var
        deltacost: real;
            { the change in the total cost of the path under the randomly 
              selection operation }
        start, finish: VERTEXRANGE;
            { the endpoints of the segment to be reversed or transposed }
        break: VERTEXRANGE;
            { the position after which a splice is to be inserted }
    begin
        case randomint (1, 2) of
            1:
                begin
                    deltacost := reversal (p, distance, start, finish);
                    changemade := Metropolis (deltacost, mutability);
                    if changemade then
                       reverse (p, start, finish, deltacost)
                end;
            2:
                begin
                    deltacost :=
                        transposal (p, distance, start, finish, break);
                    changemade := Metropolis (deltacost, mutability);
                    if changemade then
                       splice (p, start, finish, break, deltacost)
                end
        end { case }
    end;

    { writesolution

      This procedure writes a description of the completed solution
      to standard output, including its total cost }
    procedure writesolution (var problem: GRAPH; var solution: PATH;
        var distance: COSTTABLE);
    var
        city: VERTEXRANGE;
            { departure points for successive legs of the trip }
    begin
        with problem, solution do begin
            for city := 1 to vertices - 1 do begin
                write (distance[tour[city], tour[city + 1]]:10:2, '  ');
                write ('From ');
                with items[tour[city]].vname do
                    write (chars:len);
                write (' to ');
                with items[tour[city + 1]].vname do
                    write (chars:len);
                writeln
            end;
            write (distance[tour[vertices], tour[1]]:10:2, '  ');
            write ('From ');
            with items[tour[vertices]].vname do
                write (chars:len);
            write (' to ');
            with items[tour[1]].vname do
                write (chars:len);
            writeln;
            writeln ('Total cost = ', cost:10:2)
        end
    end;

begin { main program }
    randomize;
    readproblem (map);
    computedistances (map, distance);
    initialpath (map, distance, current);
    maxtrials := 500 * map.size;
    maxchanges := 10 * map.size;
    mutability := findrange (current, distance, maxtrials);
    repeat
        mutability := 0.9 * mutability;
        changes := 0;
        trials := 0;
        while (trials < maxtrials) and (changes < maxchanges) do begin
            trials := trials + 1;
            mutate (current, mutability, distance, changemade);
            if changemade then
                changes := changes + 1
        end;
        writeln ('trials = ', trials:1);
        writeln ('changes = ', changes:1);
        writeln ('mutability = ', mutability:1:4);
        writeln ('cost of current solution: ', current.cost:1:3);
        writeln
    until changes = 0;
    writesolution (map, current, distance)
end.