unit U_TSP;
 {Copyright 2002, Gary Darby, Intellitech Systems Inc., www.DelphiForFun.org

 This program may be used or modified for any non-commercial purpose
 so long as this original notice remains in place.
 All other rights are reserved
 }

 {The Traveling Salesman Problem - given a set of random US cities, find the
  shortest path connecting all of the cities.   This allows human user to draw
  their best attemnpt at a short path and then see the computer solution}
  
interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls, Spin, GIFImage, ExtCtrls;

type
  TModes=(unknown, define, path);
  TCoordObj=class(TObject)
    cityname:string;
    x,y:integer;
    visited:boolean;
  end;


  TForm1 = class(TForm)
    Image1: TImage;
    DefineBtn: TButton;
    IteneraryBtn: TButton;
    CityCountSpinEdit: TSpinEdit;
    Label1: TLabel;
    ShortestPathBtn: TButton;
    CityLbl: TLabel;
    ResetPathBtn: TButton;
    MilesLbl: TLabel;
    ShortMilesLbl: TLabel;
    RoundtripBox: TCheckBox;
    Userslist: TListBox;
    Programslist: TListBox;
    MaxtimeGrp: TRadioGroup;
    StopBtn: TButton;
    Label2: TLabel;
    Label3: TLabel;
    SearchTimeLbl: TLabel;
    Timer1: TTimer;
    Button1: TButton;
    LabelBox: TCheckBox;
    IntroBtn: TButton;
    procedure DefineBtnClick(Sender: TObject);
    procedure FormActivate(Sender: TObject);
    procedure Image1MouseUp(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure Image1MouseMove(Sender: TObject; Shift: TShiftState; X,
      Y: Integer);
    procedure IteneraryBtnClick(Sender: TObject);
    procedure ResetPathBtnClick(Sender: TObject);
    procedure ShortestPathBtnClick(Sender: TObject);
    procedure StopBtnClick(Sender: TObject);
    procedure Timer1Timer(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure CityCountSpinEditChange(Sender: TObject);
    procedure IntroBtnClick(Sender: TObject);
    procedure Button2Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    mode:TModes;
    CityList, ListToShow:TStringList;
    fname:string; {name of citylist file}
    coordobj:TCoordObj;
    oldx,oldy:integer;
    oldname:string;
    key:string;
    index:integer;
    newflag:boolean;
    bitmap:TBitmap;
    shortpathpoints, pathpoints:array of TPoint; {ordered list of cities on path}
    shortpathlength, pathlength:integer;
    scale:single;
    starttime:TDatetime; {start time of current shortest path search}
    procedure showcity(coordobj:TCoordObj);
    procedure SaveList;
    function oncity(x,y:integer; var name:string):boolean;
    procedure showcities(list:TStringlist);
    procedure drawpath;
    procedure drawshortpath;
  end;

var
  Form1: TForm1;

implementation

uses combo, U_CityDlg, U_Intro, uMakeCaption;

{$R *.DFM}

const maxpoints=51;

{*************** DefineBtnClick *************}
procedure TForm1.DefineBtnClick(Sender: TObject);
var  i:integer;
begin
  Resetpathbtnclick(sender);
   mode:=define;
   {show existing cities}
   if citylist.count>0 then
   for i:= 0 to citylist.count-1 do showcity(TCoordobj(citylist.objects[i]));
end;

{****************** FormActivate *************}
procedure TForm1.FormActivate(Sender: TObject);
var
  list:TStringList;
  i,j:integer;
  s:string;
  lapoint,nycpoint:TPoint;
begin
  makecaption('The Traveling Salesman Problem',
               #169+' 2002, G Darby, www.delphiforfun.org',self); 
  IntroForm.showmodal;
  randomize;
  mode:=unknown;
  CityList:=TstringList.Create;
  citylist.sorted:=true;
  ListToShow:=TStringList.create;
  ListToShow.sorted:=true;
  bitmap:=TBitmap.create;
  bitmap.height:=image1.height;
  bitmap.width:=image1.width;
  bitmap.loadfromfile(extractfilepath(application.exename)+'usa.bmp');
  image1.canvas.stretchdraw(rect(0,0,image1.width,image1.height),bitmap);
  fname:=extractfilepath(application.exename)+'citylist.str';
  if fileexists(fname) then
  begin
    list:=TStringlist.create;
    list.loadfromfile(fname);
    for i:=0 to list.count-1 do
    begin
      s:=list[i];
      j:=citylist.add(copy(s,1,8));

      citylist.objects[j]:=TCoordObj.create;
      with TCoordobj(citylist.objects[j]) do
      begin
        cityname:=trim(copy(s,9,32));
        x:=strtoint(copy(s,41,4));
        y:=strtoint(copy(s,45,4));
        if cityname='New York City'
        then nycpoint:=point(x,y)
        else if cityname='Los Angeles' then lapoint:=point(x,y);
      end;
      showcity(TCoordobj(citylist.objects[i]));

    end;
    list.free;

    if (nycpoint.x>0) and (lapoint.x>0) then
    scale:=3000/sqrt(sqr(nycpoint.x-lapoint.x)+sqr(nycpoint.y-lapoint.y))
    else scale:=6.7335;

  end;
  setlength(pathpoints, CityCountSpinEdit.maxvalue+2);
  setlength(shortpathpoints, CityCountSpinEdit.maxvalue+2);
  stopbtn.bringtofront; {we need the stop button in front to hide other buttons,
                         but hiding them at design time is not nice, so we do it
                         at runtime}
end;

{******************** Image1MouseMove **************}
procedure TForm1.Image1MouseMove(Sender: TObject; Shift: TShiftState; X,
  Y: Integer);
var cname:string;
begin
   if (labelbox.checked) and Oncity(x,y, cname) then
   with citylbl do
   begin
     left:=x+image1.left+5;
     top:=y+image1.top-13;
     caption:=cname;
     visible:=true;
   end
   else citylbl.visible:=false;
end;

{********************** IMage1MoueseUp ***************}
procedure TForm1.Image1MouseUp(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
   If mode=define then
   begin {define cities}
     {round coordinates to 10 pixel increments}
     key:=format('%4d%4d',[10*((x+5) div 10), 10*((y+5) div 10)]);
     if not citylist.find(key, index) then
     begin {new entry}
       coordobj:=TCoordobj.create;
       coordobj.x:=10*((x+5) div 10); {round locations}
       coordobj.y:=10*((y+5) div 10);
       coordobj.cityname:='';
       CityDlg.nameedt.text:='';
       newflag:=true;
     end
     else
     begin
       newflag:=false;
       coordobj:=tcoordobj(citylist.objects[index]);
       CityDlg.nameedt.text:=coordobj.cityname;
       oldname:=coordobj.cityname;
       oldx:=coordobj.x;
       oldy:=coordobj.y;
     end;
     CityDlg.showmodal;
   end
   else
   if mode=path then  {make a path}
   begin
     key:=format('%4d%4d',[10*((x+5) div 10), 10*((y+5) div 10)]);
     if listtoshow.find(key, index) then
     begin
       if not tCoordobj(listtoshow.objects[index]).visited then
       with tCoordobj(listtoshow.objects[index]) do
       begin
         visited:=true;
         inc(pathlength);
         pathpoints[pathlength]:=point(x,y);
         drawpath;
         userslist.items.add(cityname);
         if (roundtripbox.checked) and (pathlength=listtoshow.count)
         then {all cities visited and round trip option set - close the loop}
         begin
           visited:=true;
           inc(pathlength);
           pathpoints[pathlength]:=pathpoints[1];
           userslist.items.add(userslist.items[0]);
           drawpath;
         end;
       end;
     end;
   end;
end;

{***************** DrawPath *************}
procedure TForm1.drawpath;
var i:integer;
     miles:single;
begin

  if (mode=path) then
  begin
    miles:=0;

    if (pathlength>1) then
    with image1.canvas do
    begin
      pen.color:=clblack;

      pen.width:=3;
      pen.style:=psSolid;
      moveto(pathpoints[1].x,pathpoints[1].y);
      for i:=2 to pathlength do
      begin
        lineto(pathpoints[i].x,pathpoints[i].y);
        miles:=miles+ scale*sqrt(sqr(pathpoints[i-1].x-pathpoints[i].x)+
                           sqr(pathpoints[i-1].y-pathpoints[i].y));
      end;
    end;
    mileslbl.caption:=format('Path length %5.0f miles',[miles]);
    mileslbl.visible:=true;
  end
  else mileslbl.visible:=false;
end;

{********************* DawShortpath *************}
procedure TForm1.drawShortPath;
var i:integer;
     miles:single;
begin

  if (mode=path) then
  begin
    miles:=0;
    if (shortpathlength>1) then
    with image1.canvas do
    begin
      pen.color:=clyellow;
      brush.color:=clred;
      pen.width:=1;
      pen.style:=psdot;
      moveto(shortpathpoints[1].x,shortpathpoints[1].y);
      for i:=2 to shortpathlength do
      begin
        lineto(shortpathpoints[i].x,shortpathpoints[i].y);
        miles:=miles+ scale*sqrt(sqr(shortpathpoints[i-1].x-shortpathpoints[i].x)+
                           sqr(shortpathpoints[i-1].y-shortpathpoints[i].y));
      end;
    end;
    shortmileslbl.caption:=format('Shortest path length %5.0f miles',[miles]);
    shortmileslbl.visible:=true;
  end
  else shortmileslbl.visible:=false;
end;

{****************** ShowCity ****************}
procedure TForm1.showcity(coordobj:TCoordObj);
begin
  with coordobj, image1.picture.bitmap.canvas  do
  begin
    brush.color:=clblack;
    pen.color:=clblack;
    ellipse(x-3,y-3,x+3,y+3);
  end;
end;

{******************** ShowCities **********}
procedure TForm1.showcities(list:TStringlist);
var i:integer;
begin
  image1.canvas.stretchdraw(rect(0,0,image1.width,image1.height),bitmap);
  for i:=0 to list.count-1 do showcity(tCoordobj(list.objects[i]));
end;

{******************** SaveList **************}
procedure TForm1.saveList;
var
  i:integer;
  list:TStringlist;
  s:string;
begin
  list:=TStringlist.create;
  for i:=0 to citylist.count-1 do
  begin
    with TCoordobj(citylist.objects[i]) do
    s:=format('%8s%32s%4d%4d',[citylist[i],cityname,x,y]);
    list.add(s);
  end;
  list.savetofile(fname);
  list.free;
end;


{******************* OnCity ****************}
function TForm1.oncity(x,y:integer; var name:string):boolean;
{test if point (x,y) is on a city, if so return city name and set result true}
var
  key:string;
  index:integer;
begin
  key:=format('%4d%4d',[10*((x+5) div 10), 10*((y+5) div 10)]);
  if citylist.find(key,index) then
  begin
    name:=TCoordobj(citylist.objects[index]).cityname;
    result:=true;
  end
  else result:=false;
end;

{******************** ItineraryBtnClick ************}
procedure TForm1.IteneraryBtnClick(Sender: TObject);
var i,j,n:integer;
begin
  listtoshow.clear;
  userslist.clear;
  programslist.clear;
  i:=0;
  (*
  if citycountspinedit.value> citylist.count div 2
  then citycountspinedit.value := citylist.count div 2;
  *)
  while i< citycountspinedit.value do
  begin
    n:=random(citylist.count);
    if not listtoshow.find(citylist[n],index) then
    begin
      j:=listtoshow.addobject(citylist[n], citylist.objects[n]);
      inc(i);
      TCoordobj(listtoshow.objects[j]).visited:=false;
    end;
    showcities(listtoshow);
    mode:=path;
    pathlength:=0;
  end;
end;

{******************* ResetPathBtnClick *************}
procedure TForm1.ResetPathBtnClick(Sender: TObject);
var i:integer;
begin
  pathlength:=0;
  showcities(listtoshow);
  for i:= 0 to listtoshow.count-1 do
  tCoordobj(listtoshow.objects[i]).visited:=false;
  mileslbl.visible:=false;
  userslist.clear;
  programslist.clear;
end;

{************** ShortestPathBtnClick ************}
procedure TForm1.ShortestPathBtnClick(Sender: TObject);
var
  nbrpoints:integer;
  ShortPath: array[1..maxpoints] of integer;
  minpath:integer;
  d:array[1..maxpoints,1..maxpoints] of integer; {pairwise distances between points}
  points:array[1..maxpoints] of TPoint;


    function GetNextPath:integer;
    var
      i,dist:integer;
    Begin
      if combos.getnextpermute then
      Begin
        shortpath[1]:=combos.selected[1];
        {Initialize total path length (d) with distance from last point to first point}
        with combos do if roundtripbox.checked
        then dist:=d[selected[nbrpoints],selected[1]]
        else dist:=0;
        for i:= 2 to nbrpoints do
        Begin
          shortpath[i]:=combos.selected[i];
          with points[shortpath[i]] do
          begin
            dist:=dist+d[shortpath[i],shortpath[i-1]];
            if dist>minpath then break;  {quit if path is not shortest}
          end;
        end;
        result:=dist;
      end
      else result:=0;
    end;

    function Getfirstpath:integer;
    {Initialize shortest path length search}
    Begin
      combos.setup(nbrpoints,nbrpoints,permutations);
      minpath:=100000;
      result:=getnextpath;
    end;


var
  i,j, x1,x2,y1,y2:integer;
  pathlength:integer;
  count:integer;
  mindist, dist, lowdistindex:integer;
  temppoints: array[1..maxpoints] of TPoint;
  used: array[1..maxpoints] of boolean;
  s:string;
  first:boolean;

  maxtime:single;
begin
  If listtoshow.count<2 then exit;
  tag:=0; {reset stop flag}
  count:=0;
  case maxtimegrp.itemindex of
    0:maxtime:=15;
    1:maxtime:=30;
    2:maxtime:=60;
    3:maxtime:=300;
    else maxtime:=24*3600;
  end;

  {Rather than juggle the points for each path, we'll just keep an array
   of point numbers and use those as pointers into the points array
   to draw a path}
  screen.cursor:=crHourGlass;

  {enhancement - sort the points so that each is near its closest neighbors}
   nbrpoints:=listtoshow.count;

   with tcoordobj(listtoshow.objects[0]) do
   begin
     points[1]:=point(x,y);
     temppoints[1]:=points[1];
     used[1]:=true;
   end;
   for i:= 1 to nbrpoints-1 do
   with tcoordobj(listtoshow.objects[i]) do
   begin
     points[i+1]:=point(x,y);
     used[i+1]:=false;
   end;
   {sort the points}
   for i:= 2 to nbrpoints do
   begin
      x1:=points[i-1].x;
      y1:=points[i-1].y;
      mindist:=high(integer);
      for j:=2 to nbrpoints do
      if not used[j] then
      begin
        x2:=points[j].x;
        y2:=points[j].y;
        dist:=trunc(sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)));
        if dist<mindist then
        begin
          mindist:=dist;
          lowdistindex:=j
        end
      end;
      used[lowdistindex]:=true;
      temppoints[i]:=points[lowdistindex];
    end;
    for i:= 1 to nbrpoints do points[i]:=temppoints[i];

  {set an array of distances between each pair of points to save calc time}
  for i:= 1 to nbrpoints do
  for j:= i to nbrpoints do
  begin
    x1:=points[i].x;
    y1:=points[i].y;
    x2:=points[j].x;
    y2:=points[j].y;
    d[i,j]:=trunc(sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)));
    d[j,i]:=d[i,j];
  end;

  stopbtn.visible:=true;

  minpath:=getfirstpath;
  first:=true;
  starttime:=now;
  timer1.enabled:=true;
  repeat
    pathlength:=getnextpath; {returns 0 when no more paths}
    if ((pathlength>0) and (pathlength<minpath)) or first then
    Begin
      programslist.clear;
      minpath:=pathlength;

      for i:=1 to nbrpoints do
      begin
        shortpathpoints[i]:=points[shortpath[i]];
        with shortpathpoints[i] do if oncity(x,y, s)
        then programslist.items.add(s)
        else showmessage('City not found, tell Grandpa');
      end;
      shortpathlength:=nbrpoints;
      if roundtripbox.checked then
      begin
        inc(shortpathlength);
        shortpathpoints[shortpathlength]:=shortpathpoints[1];
        programslist.items.add(programslist.items[0]);
      end;
      showcities(listtoshow);
      drawpath;
      drawshortpath;
      application.processmessages;
      first:=false;
    end;
    inc(count);
    if count mod 8192 = 0 then
    Begin
      application.processmessages;  {handle possible stop button press}
      if tag=1 then pathlength:=0;  {stop button was hit}
      if (now-starttime)*secsperday>maxtime then pathlength:=0;
    end;
  until pathlength=0;
  screen.cursor:=crDefault;
  stopbtn.visible:=false;
  timer1.enabled:=false;
  timer1timer(self);  {make sure time display gets called at least once}
end;

{***************** StopBtnClick *************}
procedure TForm1.StopBtnClick(Sender: TObject);
begin    tag:=1;  end;

{***************** Timer1Timer *************}
procedure TForm1.Timer1Timer(Sender: TObject);
{Pops every second while exhaustive search is under way to update time display}
begin
  SearchTimeLbl.caption:=formatDateTime('"Search time (hh:mm:ss): " hh:nn:ss',
                                         now-starttime);
end;



type
  arrn=array[1..maxpoints] of integer;
  arrnn=array[1..maxpoints,1..maxpoints] of integer;


procedure TForm1.Button1Click(Sender: TObject);

    {Branch and Bounds Algorithm - not used -
               causes range check error for larger cases}
    procedure BABTSP(
           N,INF  :integer;
       var W      :ARRNN;
       var ROUTE  :ARRN;
       var TWEIGHT:integer);

       var BACKPTR,BEST,COL,FWDPTR,ROW:ARRN;
           I,INDEX                    :integer;

       procedure EXPLORE(EDGES,COST:integer;var ROW,COL:ARRN);
          var AVOID,C,COLROWVAL,FIRST,I,J,
              LAST,LOWERBOUND,MOST,R,SIZE :integer;
              COLRED,NEWCOL,NEWROW,ROWRED :ARRN;

          function MIN(I,J:integer):integer;
          begin
             if I <= J then MIN:=I else MIN :=J
          end;  { MIN }

          function REDUCE(var ROW,COL,ROWRED,COLRED:ARRN):integer;
             var I,J,RVALUE,TEMP:integer;
          begin
             RVALUE:=0;
             for I:=1 to SIZE do begin                    { REDUCE ROWS }
                TEMP:=INF;
                for J:=1 to SIZE do TEMP:=MIN(TEMP,W[ROW[I],COL[J]]);
                if TEMP > 0 then begin
                   for J:=1 to SIZE do
                      if W[ROW[I],COL[J]] < INF then
                         W[ROW[I],COL[J]]:=W[ROW[I],COL[J]]-TEMP;
                   RVALUE:=RVALUE+TEMP
                end;
                ROWRED[I]:=TEMP
             end;  { for I }
             for J:=1 to SIZE do begin                 { REDUCE COLUMNS }
                TEMP:=INF;
                for I:=1 to SIZE do TEMP:=MIN(TEMP,W[ROW[I],COL[J]]);
                if TEMP > 0 then begin
                   for I:=1 to SIZE do
                      if W[ROW[I],COL[J]] < INF then
                         W[ROW[I],COL[J]]:=W[ROW[I],COL[J]]-TEMP;
                   RVALUE:=RVALUE+TEMP
                end;
                COLRED[J]:=TEMP
             end;  { for J }
             REDUCE:=RVALUE
          end;  { REDUCE }

          procedure BESTEDGE(var R,C,MOST:integer);
             var I,J,K,MINCOLELT,MINROWELT,ZEROES:integer;
          begin
             MOST:=-INF;
             r:=0; c:=0;

             for I:=1 to SIZE do
                for J:=1 to SIZE do
                   if W[ROW[I],COL[J]] = 0 then begin
                      MINROWELT:=INF;  ZEROES:=0;
                      for K:=1 to SIZE do
                         if W[ROW[I],COL[K]] = 0 then ZEROES:=ZEROES+1
                         else MINROWELT:=MIN(MINROWELT,W[ROW[I],COL[K]]);
                      if ZEROES > 1 then MINROWELT:=0;
                      MINCOLELT:=INF;  ZEROES:=0;
                      for K:=1 to SIZE do
                         if W[ROW[K],COL[J]] = 0 then ZEROES:=ZEROES+1
                         else MINCOLELT:=MIN(MINCOLELT,W[ROW[K],COL[J]]);
                      if ZEROES > 1 then MINCOLELT:=0;
                      if (MINROWELT+MINCOLELT) > MOST then begin
                                         { A BETTER EDGE HAS BEEN FOUND }
                         MOST:=MINROWELT+MINCOLELT;
                         R:=I;  C:=J
                      end
                   end ; { if W[ROW[I],COL[J]] = 0, for J, I }
                  if (r=0) or (c=0) then showmessage('Invalid bestedge call')
          end;  { BESTEDGE }

       begin                                          { BODY OF EXPLORE }
          SIZE:=N-EDGES;
          COST:=COST+REDUCE(ROW,COL,ROWRED,COLRED);
          if COST < TWEIGHT then
             if EDGES = (N-2) then
             begin    { LAST TWO EDGES ARE FORCED }
                for I:=1 to N do BEST[I]:=FWDPTR[I];
                if W[ROW[1],COL[1]] = INF then AVOID:=1
                else AVOID:=2;
                BEST[ROW[1]]:=COL[3-AVOID];  BEST[ROW[2]]:=COL[AVOID];
                TWEIGHT:=COST
             end  { if EDGES = (N-2) }
             else
             begin
                BESTEDGE(R,C,MOST);
                LOWERBOUND:=COST+MOST;
                FWDPTR[ROW[R]]:=COL[C];  BACKPTR[COL[C]]:=ROW[R];
                LAST:=COL[C];       { PREVENT CYCLES }
                while FWDPTR[LAST] <> 0 do LAST:=FWDPTR[LAST];
                FIRST:=ROW[R];
                while BACKPTR[FIRST] <> 0 do FIRST:=BACKPTR[FIRST];
                COLROWVAL:=W[LAST,FIRST];  W[LAST,FIRST]:=INF;
                for I:=1 to R-1 do NEWROW[I]:=ROW[I];      { REMOVE ROW }
                for I:=R to SIZE-1 do NEWROW[I]:=ROW[I+1];
                for I:=1 to C-1 do NEWCOL[I]:=COL[I];      { REMOVE COL }
                for I:=C to SIZE-1 do NEWCOL[I]:=COL[I+1];
                EXPLORE(EDGES+1,COST,NEWROW,NEWCOL);
                W[LAST,FIRST]:=COLROWVAL;     { RESTORE PREVIOUS VALUES }
                BACKPTR[COL[C]]:=0;  FWDPTR[ROW[R]]:=0;
                if LOWERBOUND < TWEIGHT then
                begin
                   W[ROW[R],COL[C]]:=INF;
                   EXPLORE(EDGES,COST,ROW,COL);
                   W[ROW[R],COL[C]]:=0
                end
             end;  { else: EDGES < N-2 }
          for I:=1 to SIZE do                         { UNREDUCE MATRIX }
             for J:=1 to SIZE do
                W[ROW[I],COL[J]]:=W[ROW[I],COL[J]]+ROWRED[I]+COLRED[J]
       end;  { EXPLORE }

    begin                                                   { MAIN BODY }
       for I:=1 to N do
       begin
          ROW[I]:=I;  COL[I]:=I;
          FWDPTR[I]:=0;  BACKPTR[I]:=0
       end;
       TWEIGHT:=INF;
       EXPLORE(0,0,ROW,COL);
       INDEX:=1;
       for I:=1 to N do
       begin
          ROUTE[I]:=INDEX;  INDEX:=BEST[INDEX]
       end
    end;  { BABTSP }

  {**********************************************}
  {Furthest Insertion TSP}
   PROCEDURE FITSP(N,S,INF:INTEGER; VAR W:ARRNN; VAR ROUTE:ARRN; VAR TWEIGHT:INTEGER);

   VAR END1,END2,FARTHEST,I,INDEX,
          INSCOST,J,MAXDIST,NEWCOST,NEXTINDEX :INTEGER;
       CYCLE,DIST                             :ARRN;
BEGIN
   FOR I:=1 TO N DO CYCLE[I]:=0;
   CYCLE[S]:=S;
   FOR I:=1 TO N DO DIST[I]:=W[S,I];
   TWEIGHT:=0;
   FOR I:=1 TO N-1 DO
   BEGIN
      MAXDIST:=-INF;
      FOR J:=1 TO N DO
         IF CYCLE[J] = 0 THEN
            IF DIST[J] > MAXDIST THEN
            BEGIN
               MAXDIST:=DIST[J];  FARTHEST:=J
            END;
      INSCOST:=INF;  INDEX:=S;
      FOR J:=1 TO I DO
      BEGIN
         NEXTINDEX:=CYCLE[INDEX];
         NEWCOST:=W[INDEX,FARTHEST]+W[FARTHEST,NEXTINDEX]-
                  W[INDEX,NEXTINDEX];
         IF NEWCOST < INSCOST THEN
         BEGIN
            INSCOST:=NEWCOST;
            END1:=INDEX;  END2:=NEXTINDEX
         END;
         INDEX:=NEXTINDEX
      END;  { FOR J }
      CYCLE[FARTHEST]:=END2;  CYCLE[END1]:=FARTHEST;
      TWEIGHT:=TWEIGHT+INSCOST;
      FOR J:=1 TO N DO
         IF CYCLE[J] = 0 THEN
            IF W[FARTHEST,J] < DIST[J] THEN DIST[J]:=W[FARTHEST,J]
   END;  { FOR I }
   INDEX:=S;
   FOR I:=1 TO N DO
   BEGIN
      ROUTE[I]:=INDEX;  INDEX:=CYCLE[INDEX]
   END
END;  { FITSP }

{****************************************************}
{Two point exchange optimization}
PROCEDURE TWOOPT( N:INTEGER; VAR W:ARRNN; VAR ROUTE:ARRN; VAR TWEIGHT:INTEGER);

   VAR AHEAD,I,I1,I2,INDEX,J,J1,J2,
       LAST,LIMIT,MAX,MAX1,NEXT,S1,S2,T1,T2:INTEGER;
       PTR                                 :ARRN;
BEGIN
   FOR I:=1 TO N-1 DO PTR[ROUTE[I]]:=ROUTE[I+1];
   PTR[ROUTE[N]]:=ROUTE[1];
   REPEAT  { UNTIL MAX = 0 }
      MAX:=0;  I1:=1;
      FOR I:=1 TO N-2 DO
      BEGIN
         IF I = 1 THEN LIMIT:=N-1
         ELSE LIMIT:=N;
         I2:=PTR[I1];  J1:=PTR[I2];
         FOR J:=I+2 TO LIMIT DO
         BEGIN
            J2:=PTR[J1];
            MAX1:=W[I1,I2]+W[J1,J2]-(W[I1,J1]+W[I2,J2]);
            IF MAX1 > MAX THEN
            BEGIN   { BETTER PAIR HAS BEEN FOUND }
               S1:=I1;  S2:=I2;
               T1:=J1;  T2:=J2;
               MAX:=MAX1
            END;
            J1:=J2
         END;  { FOR J }
         I1:=I2
      END;  { FOR I }
      IF MAX > 0 THEN
      BEGIN                    { SWAP PAIR OF EDGES }
         PTR[S1]:=T1;
         NEXT:=S2;  LAST:=T2;
         REPEAT
            AHEAD:=PTR[NEXT];  PTR[NEXT]:=LAST;
            LAST:=NEXT;  NEXT:=AHEAD
         UNTIL NEXT = T2;
         TWEIGHT:=TWEIGHT-MAX
      END  { IF MAX > 0 }
   UNTIL MAX = 0;
   INDEX:=1;
   FOR I:=1 TO N DO
   BEGIN
      ROUTE[I]:=INDEX;  INDEX:=PTR[INDEX]
   END
END;  { TWOOPT }

{***************************************************}


{three point exchange optimization}
procedure THREEOPT(N:integer; var W:ARRNN; var ROUTE:ARRN);

   type SWAPTYPE  =(ASYMMETRIC,SYMMETRIC);
        SWAPRECORD=record
                      X1,X2,Y1,Y2,Z1,Z2,GAIN:integer;
                      CHOICE                :SWAPTYPE
                   end;
   var BESTSWAP,SWAP:SWAPRECORD;
       I,INDEX,J,K  :integer;
       PTR          :ARRN;

   procedure SWAPCHECK(var SWAP:SWAPRECORD);
      var DELWEIGHT,MAX:integer;
   begin
      with SWAP do begin
         GAIN:=0;
         DELWEIGHT:=W[X1,X2]+W[Y1,Y2]+W[Z1,Z2];
         MAX:=DELWEIGHT-(W[Y1,X1]+W[Z1,X2]+W[Z2,Y2]);
         if MAX > GAIN then begin
            GAIN:=MAX;  CHOICE:=ASYMMETRIC
         end;
         MAX:=DELWEIGHT-(W[X1,Y2]+W[Z1,X2]+W[Y1,Z2]);
         if MAX > GAIN then begin
            GAIN:=MAX;  CHOICE :=SYMMETRIC
         end
      end  { with SWAP }
   end;  { SWAPCHECK }

   procedure REVERSE(START,FINISH:integer);
      var AHEAD,LAST,NEXT:integer;
   begin
      if START <> FINISH then begin
         LAST:=START;  NEXT:=PTR[LAST];
         repeat
            AHEAD:=PTR[NEXT];  PTR[NEXT]:=LAST;
            LAST:=NEXT;  NEXT:=AHEAD;
         until LAST = FINISH
      end  { if START <> FINISH }
   end;  { RESERVE }

begin                                                   { MAIN BODY }
   for I:=1 to N-1 do PTR[ROUTE[I]]:=ROUTE[I+1];
   PTR[ROUTE[N]]:=ROUTE[1];
   repeat  { until BESTSWAP.GAIN = 0 }
      BESTSWAP.GAIN:=0;
      SWAP.X1:=1;
      for I:=1 to N do begin
         SWAP.X2:=PTR[SWAP.X1];  SWAP.Y1:=SWAP.X2;
         for J:=2 to N-3 do begin
            SWAP.Y2:=PTR[SWAP.Y1];  SWAP.Z1:=PTR[SWAP.Y2];
            for K:=J+2 to N-1 do begin
               SWAP.Z2:=PTR[SWAP.Z1];
               SWAPCHECK(SWAP);
               if SWAP.GAIN > BESTSWAP.GAIN then BESTSWAP:=SWAP;
               SWAP.Z1:=SWAP.Z2
            end;  { for K }
            SWAP.Y1:=SWAP.Y2
         end;  { for J }
         SWAP.X1:=SWAP.X2
      end;  { for I }
      if BESTSWAP.GAIN > 0 then
         with BESTSWAP do begin
            if CHOICE = ASYMMETRIC then begin
               REVERSE(Z2,X1);
               PTR[Y1]:=X1;  PTR[Z2]:=Y2
            end
            else begin
               PTR[X1]:=Y2;  PTR[Y1]:=Z2
            end;
            PTR[Z1]:=X2;
         end  { with BESTSWAP }
   until BESTSWAP.GAIN = 0;
   INDEX:=1;
   for I:=1 to N do begin
      ROUTE[I]:=INDEX;  INDEX:=PTR[INDEX]
   end
end;  { THREEOPT }


 var
  nbrpoints:integer;
  w:arrNN; {pairwise distances between points}
  points:array[1..maxpoints] of TPoint;
  route:arrn;

  i,j, x1,x2,y1,y2:integer;
  pathlength:integer;
  s:string;
begin
  If listtoshow.count<2 then exit;
  tag:=0; {reset stop flag}
  screen.cursor:=crHourGlass;
  nbrpoints:=listtoshow.count;

   for i:= 0 to nbrpoints-1 do
   with tcoordobj(listtoshow.objects[i]) do
   begin
     points[i+1]:=point(x,y);
     route[i+1]:=i+1;
   end;
   {set an array of distances between each pair of points to save calc time}
   for i:= 1 to nbrpoints do
   for j:= i to nbrpoints do
   begin
     x1:=points[i].x;
     y1:=points[i].y;
     x2:=points[j].x;
     y2:=points[j].y;
     w[i,j]:=trunc(sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)));
     if i=j then w[i,j]:=0 {999999 for BABTSP}
     else w[j,i]:=w[i,j];
   end;

   stopbtn.visible:=true;
   starttime:=now;
   timer1.enabled:=true;

  try
    {BABTSP( Nbrpoints, maxint, W, route, pathlength);}
    FITSP(Nbrpoints,1,999999,W,ROUTE,pathlength);
    TWOOPT(nbrpoints,W,ROUTE,pathlength);
    THREEOPT(nbrpoints,W,ROUTE);
  finally
    programslist.clear;
    for i:=1 to nbrpoints do shortpathpoints[i]:=points[route[i]];
    for i:=1 to nbrpoints do
    begin
      with shortpathpoints[i] do if oncity(x,y, s)
      then programslist.items.add(s)
      else showmessage('City not found, tell Grandpa');
    end;
    shortpathlength:=nbrpoints;
    if roundtripbox.checked then
    begin
      inc(shortpathlength);
      shortpathpoints[shortpathlength]:=shortpathpoints[1];
      programslist.items.add(programslist.items[0]);
    end;
    showcities(listtoshow);
    drawpath;
    drawshortpath;
    screen.cursor:=crDefault;
    stopbtn.visible:=false;
    timer1.enabled:=false;
    timer1timer(self);  {make sure time display gets called at least once}
  end;
end;




procedure TForm1.CityCountSpinEditChange(Sender: TObject);
begin  resetpathbtnclick(sender); end;

procedure TForm1.IntroBtnClick(Sender: TObject);
begin    Introform.showmodal; end;

procedure TForm1.Button2Click(Sender: TObject);
begin
  image1.picture.Bitmap.pixelformat:=pf24bit;
  image1.picture.savetofile('TSP.bmp');
end;

end.