Programmer 7500/MAX_PROGRAMMERS.iso/CLIPPER/MISC/NETSTUFF.ZIP/TSP.PAS

Pascal/Delphi Source File | 1980-01-04 | 10.4 KB | 348 lines

{$R+}
PROGRAM traveling_salesperson ;

(* Copyright 1987 - Knowledge Garden Inc.
                    473A Malden Bridge Rd.
                    R.D. 2
                    Nassau, NY 12123       *)


(* TSP solves a series of differential equations which simulate a neural
   net solution of the traveling salesperson problem. The problem and
   the equations are described in the article "Computing with Neurons" in
   the July 1987 issue of AI Expert Magazine.

   This program has been tested using Turbo ver 3.01A on an IBM PC/AT. It has
   been run under both DOS 3.2 and Concurrent 5.0 .

   We would be pleased to hear your comments, good or bad, or any applications
   and modifications of the program. Contact us at:

     AI Expert
     500 Howard St.
     San Francisco, CA 94105

   Bill and Bev Thompson    *)

 CONST
  max_city = 'E' ;         (* max_city and max_position are the size of the *)
  max_position = 5 ;       (* neural net. They must match. Cities run from *)
                           (* A to max_city *)

  a = 500.0 ;              (* these are the weighting constants described *)
  b = 500.0 ;              (* in the article. By changing then you can *)
  c = 200.0 ;              (* get different types of solutions *)
  d = 300.0 ;              (* d seems to have the most effect, increasing *)
                           (* it produces shorter distance routes, but *)
                           (* they aren't necessarily real tours. *)

  u0 = 0.02 ;              (* This parameter effects the output voltage of *)
                           (* the amplifiers. Increasing it gives a broader *)
                           (* curve. *)

  n = 7 ;                  (* This term affects global inhibition of the *)
                           (* network. By setting it slightly larger than *)
                           (* the number of cities, we seem to get better *)
                           (* results *)

  h = 0.01 ;               (* The time step *)

 TYPE
  cities = 'A' .. max_city ;
  positions = 1 .. max_position ;


 VAR
  u : ARRAY [cities,positions] OF real ;      (* Input voltages *)
  dist : ARRAY [cities,cities] OF real ;      (* Distances between cities *)



 FUNCTION v(city : cities ; position : positions) : real ;
  (* This function calculates the output voltage from an amplifier
     tanh calculates the hyperbolic tangent which gives the shape
     of the output curve described in the article *)

  FUNCTION tanh(r : real) : real ;
   VAR
    r1,r2 : real ;
   BEGIN
    IF r > 20.0
     THEN tanh := 1.0
    ELSE IF r < -20.0
     THEN tanh := -1.0
    ELSE
     BEGIN
      r1 := exp(r) ;
      r2 := exp(-r) ;
      tanh := (r1 - r2) / (r1 + r2) ;
     END ;
   END ; (* tanh *)

  BEGIN
   v := (1.0 + tanh(u[city,position] / u0)) / 2.0 ;
  END ; (* v *)


 FUNCTION f(city : cities ; position : positions) : real ;
  (* This function calculates the right hand side of the differential
     equations described in the article. It is not optimized for anything
     and is pretty slow. *)

  FUNCTION col_sum(cty : cities) : real ;
   (* column inhibition. This function helps keep the number of
      output items in each column small *)
   VAR
    col : positions ;
    sum : real ;
   BEGIN
    sum := 0.0 ;
    FOR col := 1 TO max_position DO
     IF col <> position
      THEN sum := sum + v(cty,col) ;
    col_sum := sum ;
   END ; (* col_sum *)
è
  FUNCTION row_sum(p : positions) : real ;
   (* row inhibition. This function helps keep the number of
      output items in each row small *)
   VAR
    row : cities ;
    sum : real ;
   BEGIN
    sum := 0.0 ;
    FOR row := 'A' TO max_city DO
     IF row <> city
      THEN sum := sum + v(row,p) ;
    row_sum := sum ;
   END ; (* row_sum *)

  FUNCTION matrix_sum : real ;
   (* global inhibition. This function keeps the total number of cities
      visited small *)
   VAR
    row : cities ;
    col : positions ;
    sum : real ;
   BEGIN
    sum := 0.0 ;
    FOR row := 'A' TO max_city DO
     FOR col := 1 TO max_position DO
      sum := sum + v(row,col) ;
    matrix_sum := sum ;
   END ; (* matrix_sum *)

  FUNCTION dist_sum : real ;
   (* distance inhibition. The inhibition is larger for longer tours.
      Note that neuron (X,max_position) is connected to neuron (X,1),
      in other words, the net is circular *)
   VAR
    c : cities ;
    sum : real ;
   BEGIN
    sum := 0.0 ;
    IF position = max_position
     THEN
      FOR c := 'A' TO max_city DO
       sum := sum + dist[city,c] * (v(c,1) + v(c,position - 1))
    ELSE IF position = 1
     THEN
      FOR c := 'A' TO max_city DO
       sum := sum + dist[city,c] * (v(c,position + 1) + v(c,max_position))
    ELSE
     FOR c := 'A' TO max_city DO
      sum := sum + dist[city,c] * (v(c,position + 1) + v(c,position - 1)) ;
    dist_sum := sum ;
   END ; (* dist_sum *)

  BEGIN
   f := -u[city,position] - a * col_sum(city) - b * row_sum(position)
è        - c * (matrix_sum - n) - d * dist_sum ;
  END ; (* f *)


 PROCEDURE iterate ;
  (* The basic solution process. This is a terrible way to solve differential
     equations. Don't use it for anything serious, it performs poorly
     when the number of cities gets larger than 7 or 8.
     We keep iterating until the norm is less than tol or until the user
     gets bored and presses the space bar. *)
  CONST
   tol = 1.0E-05 ;
  VAR
   step : integer ;
   c1 : cities ;
   i : positions ;
   nr : real ;
   u_old : ARRAY [cities,positions] OF real ;
   ch : char ;

  FUNCTION norm : real ;
   (* The norm is a measure of how much change there has been between
      solutions. This is an infinity norm, calculated as the maximum
      absolute value of the difference between components of the
      solution vectors. We calculate the relative norm as:
        N(u_new - u) / N(u). *)
   VAR
    cx : cities ;
    ix : positions ;
    max,max_comp : real ;
   BEGIN
    max := 0.0 ;
    FOR cx := 'A' TO max_city DO
     FOR ix := 1 TO max_position DO
      BEGIN
       IF abs(u_old[cx,ix] - u[cx,ix]) > max
        THEN max := abs(u_old[cx,ix] - u[cx,ix]) ;
       IF abs(u[cx,ix]) > max_comp
        THEN max_comp := abs(u[cx,ix]) ;
      END ;
    norm := max / max_comp ;
   END ; (* norm *)

  PROCEDURE print_matrix ;
   (* Every so often, we print the input and output matrices so that
      you can see what is going on. If the output matrix describes a
      valid tour, we print that also. *)
   VAR
    c1 : cities ;
    i : positions ;
    vv : real ;
    t : ARRAY [1 .. max_position] OF char ;
    t_count : integer ;

   PROCEDURE write_tour ;
è    VAR
     i : positions ;
     t_dist : real ;
    BEGIN
     t_dist := 0.0 ;
     FOR i := 1 TO max_position - 1 DO
      t_dist := t_dist + dist[t[i],t[i+1]] ;
     t_dist := t_dist + dist[t[max_position],t[1]] ;
     write(output,'Tour: ') ;
     FOR i := 1 TO max_position DO
      write(output,t[i]) ;
     writeln(output,'   dist = ',t_dist) ;
    END ; (* write_tour *)

   PROCEDURE matrix_heading ;
    VAR
     i : positions ;
    BEGIN
     write(output,'  ') ;
     FOR i := 1 TO max_position DO
      write(output,i : 12) ;
     writeln ;
    END ; (* matrix_heading *)

   BEGIN
    t_count := 0 ;
    FOR i := 1 TO max_position DO
     t[i] := chr(0) ;
    writeln(output) ;
    writeln(output,'Step: ',step,' norm = ',nr) ;
    writeln(output) ;
    writeln(output,'Input Voltages') ;
    matrix_heading ;
    FOR c1 := 'A' TO max_city DO
     BEGIN
      write(output,c1,'    ') ;
      FOR i := 1 TO max_position DO
       write(output,u[c1,i] : 12 : 5) ;
      writeln(output) ;
     END ;
    writeln(output) ;
    writeln(output,'Output Voltages') ;
    matrix_heading ;
    FOR c1 := 'A' TO max_city DO
     BEGIN
      write(output,c1,'    ') ;
      FOR i := 1 TO max_position DO
       BEGIN
        vv := v(c1,i) ;
        write(output,vv : 12 : 5) ;
        IF (vv > 0.8) AND (t_count < max_position) AND (t[i] = chr(0))
         THEN
          BEGIN
           t_count := t_count + 1 ;
           t[i] := c1 ;
è          END ;
       END ;
      writeln(output) ;
     END ;
    IF t_count = max_position
     THEN write_tour ;
   END ; (* print_matrix *)

  BEGIN
   step := 0 ;
   REPEAT
    step := step + 1 ;
    move(u,u_old,sizeof(u)) ;
    FOR c1 := 'A' TO max_city DO
     FOR i := 1 TO max_position DO
      u[c1,i] := u[c1,i] + h * f(c1,i) ;
    nr := norm ;
    IF ((step MOD 10) = 0) OR (step < 10)
     THEN print_matrix ;
   UNTIL keypressed OR (nr < tol) ;
   IF keypressed
    THEN read(kbd,ch) ;
   print_matrix ;
  END ; (* iterate *)


 PROCEDURE initialize ;
  TYPE
   location = RECORD
               x : real ;
               y : real ;
              END ;
   city_array = ARRAY [cities] OF location ;
  CONST
   u00 = -0.01386 ;
(* city_loc : city_array = ( (x : 0.21192 ; y : 0.54866),
                             (x : 0.98817 ; y : 0.68465),
                             (x : 0.53109 ; y : 0.72173),
                             (x : 0.31459 ; y : 0.79397),
                             (x : 0.63290 ; y : 0.85573)) ;

   These are the values we used for the article, if you want to
   check our results, remove the comments here and use this data *)
  VAR
   c1,c2 : cities ;
   i : positions ;
   city_loc : city_array ;
   ch : char ;
  BEGIN
   randomize ;
   FOR c1 := 'A' TO max_city DO
    BEGIN
     city_loc[c1].x := random ;
     city_loc[c1].y := random ;
    END ;
è   FOR c1 := 'A' TO pred(max_city) DO
    BEGIN
     dist[c1,c1] := 0.0 ;
     FOR c2 := succ(c1) TO max_city DO
      BEGIN
       dist[c1,c2] := sqrt(sqr(city_loc[c1].x - city_loc[c2].x) +
                           sqr(city_loc[c1].y - city_loc[c2].y)) ;
       dist[c2,c1] := dist[c1,c2] ;
      END ;
    END ;
   dist[max_city,max_city] := 0.0 ;
   FOR c1 := 'A' TO max_city DO
    FOR i := 1 TO max_position DO
     u[c1,i] := u00 + (((2 * random - 1.0) / 10.0) * u0) ;
   clrscr ;
   writeln('TSP         [c] 1987 Knowledge Garden Inc.') ;
   writeln('                     473A Malden Bridge Rd') ;
   writeln('                     Nassau, NY 12123') ;
   writeln ;
   writeln('Press <Space Bar> to begin - Press again to stop iterating.') ;
   read(kbd,ch) ;
  END ; (* initialize *)


 BEGIN
  initialize ;
  iterate ;
 END.