(*^ ::[ frontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.1"; macintoshStandardFontEncoding; paletteColors = 128; alwaysUseMacCustomDithers; currentKernel; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, R65535, e8, 28, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, cellOutline, groupLikeInput, M42, N23, bold, L-5, r58981, g58981, b58981, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; ; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, B65535, 10, "Geneva"; ; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = leftheader, inactive, L2, 12, "Times"; ; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; ; fontset = leftfooter, inactive, L2, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 9, "Courier"; ; ] :[font = title; inactive; Cclosed; preserveAspect; startGroup; ] Solving the Traveling Salesman Problem using InterCall :[font = subsubtitle; inactive; preserveAspect; endGroup; ] Terry Robb :[font = section; inactive; Cclosed; preserveAspect; startGroup; ] Abstract :[font = text; inactive; preserveAspect; ] The traveling salesman problem is a classic example of combinatorial minimization. It involves finding the shortest itinerary for a traveling salesman who must visit each of N cities in turn. One technique for solving this problem is the method of simulated annealing, and this notebook uses InterCall to access an external subroutine which implements the annealing algorithm. An example in which an additional constraint, such as adding a penalty for a river crossing, will also be demonstrated. :[font = text; inactive; preserveAspect; ] For a concrete example we use the following collection of 'cities': :[font = postscript; PostScript; formatAsPostScript; output; inactive; preserveAspect; center; pictureLeft = 98; pictureWidth = 282; pictureHeight = 282; endGroup; ] %! %%Creator: Mathematica %%AspectRatio: 1 MathPictureStart %% Graphics /Courier findfont 10 scalefont setfont % Scaling calculations 0.0238095 0.960061 0.0238095 0.953334 [ [ 0 0 0 0 ] [ 1 1 0 0 ] ] MathScale % Start of Graphics 1 setlinecap 1 setlinejoin newpath [ ] 0 setdash 0 g p P 0 0 m 1 0 L 1 1 L 0 1 L closepath clip newpath p p p .03 w .52976 .87418 Mdot .38959 .65873 Mdot .46544 .449 Mdot .25902 .36796 Mdot .21582 .90183 Mdot .45104 .49571 Mdot .68337 .4509 Mdot .83026 .20685 Mdot .46736 .55291 Mdot .1515 .97619 Mdot .97619 .83414 Mdot .52592 .39084 Mdot .09197 .93138 Mdot .4376 .92662 Mdot .58833 .54814 Mdot .35599 .78648 Mdot .96275 .85512 Mdot .26478 .4223 Mdot .54032 .95617 Mdot .34831 .62536 Mdot .63921 .1058 Mdot .26575 .64348 Mdot .33967 .82938 Mdot .02573 .90564 Mdot .1707 .54719 Mdot .76498 .73214 Mdot .70545 .02762 Mdot .28879 .14202 Mdot .22254 .03144 Mdot .39055 .86751 Mdot .11982 .57674 Mdot .05741 .51096 Mdot .4232 .40324 Mdot .07085 .45281 Mdot .21582 .74739 Mdot .90227 .60153 Mdot .78898 .77504 Mdot .84754 .56244 Mdot .75154 .1201 Mdot .91571 .33555 Mdot .60849 .93234 Mdot .85618 .59676 Mdot .37135 .30123 Mdot .92915 .20208 Mdot .30607 .87609 Mdot .58257 .18397 Mdot .10925 .11056 Mdot .14478 .81603 Mdot .35311 .34413 Mdot .32143 .4347 Mdot P P P % End of Graphics MathPictureEnd :[font = section; inactive; Cclosed; pageBreak; preserveAspect; startGroup; animationSpeed = 0; ] The problem :[font = text; inactive; preserveAspect; animationSpeed = 0; ] The "traveling salesman problem" involves finding the shortest itinerary for a traveling salesman who must visit each of N cities in turn. For a concrete example we shall look at a 50 city problem. :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Setup the position of the cities :[font = text; inactive; preserveAspect; animationSpeed = 0; ] Here are 50 random positions of some cities. :[font = input; preserveAspect; animationSpeed = 0; ] positions = { {.527,.892}, {.381,.666}, {.460,.446}, {.245,.361}, {.200,.921}, {.445,.495}, {.687,.448}, {.840,.192}, {.462,.555}, {.133,.999}, {.992,.850}, {.523,.385}, {.071,.952}, {.431,.947}, {.588,.550}, {.346,.800}, {.978,.872}, {.251,.418}, {.538,.978}, {.338,.631}, {.641,.086}, {.252,.650}, {.329,.845}, {.002,.925}, {.153,.549}, {.772,.743}, {.710,.004}, {.276,.124}, {.207,.008}, {.382,.885}, {.100,.580}, {.035,.511}, {.416,.398}, {.049,.450}, {.200,.759}, {.915,.606}, {.797,.788}, {.858,.565}, {.758,.101}, {.929,.327}, {.609,.953}, {.867,.601}, {.362,.291}, {.943,.187}, {.294,.894}, {.582,.168}, {.089,.091}, {.126,.831}, {.343,.336}, {.310,.431}}; {x,y} = Transpose[positions]; :[font = text; inactive; preserveAspect; animationSpeed = 0; ] Graphically the positions of the cities, given in the abstract earlier, could be displayed with :[font = input; preserveAspect; animationSpeed = 0; ] (* cities= {PointSize[0.03], Point/@positions}; Show[ Graphics[{cities}], AspectRatio->1 ]; *) :[font = text; inactive; preserveAspect; endGroup; endGroup; ] We could use ListPlot too, but for consistency with later examples we use more basic graphic primitives instead. ;[s] 3:0,0;13,1;21,0;113,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,10,Courier,1,12,0,0,0; :[font = section; inactive; Cclosed; pageBreak; preserveAspect; startGroup; animationSpeed = 0; ] The method :[font = text; inactive; preserveAspect; animationSpeed = 0; ] We shall use the method of simulated annealing with an annealing schedule as described in the Numerical Recipes book (see references at end). Annealing is very CPU intensive (although it is one of the best methods available) and would be inefficient to code in Mathematica directly. Luckily the Numerical Recipes book gives a fortran implementation that we can import into Mathematica fairly easily using InterCall. ;[s] 5:0,0;261,1;272,0;373,1;384,0;416,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,2,12,0,0,0; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Startup InterCall and import an ANNEAL routine :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Load in the InterCall package :[font = text; inactive; preserveAspect; animationSpeed = 0; ] First we must load in the following package. :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] <In, $Y->In, $IORDER->Range[$NCITY] :> Out, $NCITY->Min[ROWS[$X],ROWS[$Y]], $ALEN->Function[{x1,x2,y1,y2}, Sqrt[(x2-x1)^2+(y2-y1)^2]], $T->0.5, $TFACTR->0.9, $NOVER->100*$NCITY, $NLIMIT->10*$NCITY, $NTSTEPS->100, $MONITR->False, $IWRK:>Null ], S[R[$NCITY],R[$NCITY],I[$NCITY],I, RF[R,R,R,R],R,R,I,I,I,L,I[$NCITY]], "./anneal.o" ]; :[font = text; inactive; preserveAspect; animationSpeed = 0; ] You should read the Numerical Recipes book to understand the arguments and default settings indicated above. The default setting says that the default cost function ALEN between any two cities is the euclidean distance between those cities. :[font = text; inactive; preserveAspect; endGroup; ] As in the book the annealing schedule is based on the following: T is decreased by TFACTR for at most NTSTEPS, and at each step it is held constant for NOVER reconfigurations or NLIMIT successful reconfigurations whichever comes first. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Check the default setting :[font = text; inactive; preserveAspect; animationSpeed = 0; ] Having entered a default setting for the anneal routine, we can inspect the setting by :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] GetDefault[anneal] :[font = print; inactive; preserveAspect; animationSpeed = 0; ] ANNEAL[ (* TYPE=S *) $X -> In, (* DATA=R[$NCITY] *) $Y -> In, (* DATA=R[$NCITY] *) $IORDER -> Range[$NCITY] :> Out, (* DATA=I[$NCITY] *) $NCITY -> Min[ROWS[$X], ROWS[$Y]], (* DATA=I *) $ALEN -> Function[{x1, x2, y1, y2}, Sqrt[(x2 - x1)^2 + (y2 - y1)^2]], (* DATA=RF[R, R, R, R] *) $T -> 0.5, (* DATA=R *) $TFACTR -> 0.9, (* DATA=R *) $NOVER -> 100*$NCITY, (* DATA=I *) $NLIMIT -> 10*$NCITY, (* DATA=I *) $NTSTEPS -> 100, (* DATA=I *) $MONITR -> False, (* DATA=L *) $IWRK :> Null (* DATA=I[$NCITY] *) ] (* CODE="./anneal.o" *) ;[s] 1:0,1;803,-1; 2:0,12,10,Courier,0,12,0,0,0;1,11,9,Courier,0,10,0,0,0; :[font = output; output; inactive; preserveAspect; endGroup; animationSpeed = 0; ] DefaultEntry[anneal, ANNEAL, {{$X, Hold[In], None}, {$Y, Hold[In], None}, {$IORDER, Hold[Range[$NCITY]], Keep[Out, #1 & ]}, {$NCITY, Hold[Min[ROWS[$X], ROWS[$Y]]], None}, {$ALEN, Hold[Function[{x1, x2, y1, y2}, Sqrt[(x2 - x1)^2 + (y2 - y1)^2]]], None}, {$T, Hold[0.5], None}, {$TFACTR, Hold[0.9], None}, {$NOVER, Hold[100*$NCITY], None}, {$NLIMIT, Hold[10*$NCITY], None}, {$NTSTEPS, Hold[100], None}, {$MONITR, Hold[False], None}, {$IWRK, None, Keep[Null, #1 & ]}}, S, {R[$NCITY], R[$NCITY], I[$NCITY], I, RF[R, R, R, R], R, R, I, I, I, L, I[$NCITY]}, "../anneal.o", {1, 2, 4}] ;[o] anneal[$X_, $Y_] -> $IORDER :[font = text; inactive; preserveAspect; endGroup; animationSpeed = 0; ] The summary line at the end of the output above, shows that we need only supply the co-ordinates X and Y of the cities, and that the best order IORDER will be returned. (The summary line alone could be obtained by typing just ?anneal.) ;[s] 3:0,0;226,1;233,0;236,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,10,Courier,1,12,0,0,0; :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Import the routine into Mathematica ;[s] 3:0,0;24,1;35,0;37,-1; 2:2,13,9,Times,1,12,0,0,0;1,13,9,Times,3,12,0,0,0; :[font = text; inactive; preserveAspect; animationSpeed = 0; ] To import the routine, and any other routines for that matter, into Mathematica we just use InterCall's Import command. ;[s] 5:0,0;68,1;79,0;104,2;110,0;120,-1; 3:3,13,9,Times,0,12,0,0,0;1,13,9,Times,2,12,0,0,0;1,13,10,Courier,1,12,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Import[{ANNEAL}]; :[font = message; inactive; preserveAspect; fontSize = 10; animationSpeed = 0; ] InterCall::opened: Opened connection to host solar :[font = message; inactive; preserveAspect; fontSize = 10; animationSpeed = 0; ] InterCall::import: Importing: {ANNEAL} :[font = message; inactive; preserveAspect; fontSize = 10; endGroup; animationSpeed = 0; ] InterCall::linked: Using remote driver version 2.0 on host solar :[font = text; inactive; preserveAspect; endGroup; endGroup; animationSpeed = 0; ] In our case the output indicates that a connection is made to the computer named solar which in this case is the machine running the Mathematica kernel. It's possible to import code from another machine too however. ;[s] 3:0,0;133,1;144,0;216,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Call ANNEAL to find the best order to visit all the cities :[font = text; inactive; preserveAspect; animationSpeed = 0; ] For this example, the ALEN function is called about a million times, and so for efficiency the InterCall compiler/interpreter must be switched on. To compute a near optimal path the following takes about 30 seconds of real time to execute. :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] order = anneal[x,y] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {36, 42, 38, 40, 44, 8, 39, 27, 21, 46, 12, 7, 15, 9, 6, 3, 33, 50, 18, 4, 49, 43, 28, 29, 47, 34, 32, 31, 25, 22, 20, 2, 35, 48, 24, 13, 10, 5, 45, 23, 16, 30, 14, 1, 19, 41, 26, 37, 17, 11} ;[o] {36, 42, 38, 40, 44, 8, 39, 27, 21, 46, 12, 7, 15, 9, 6, 3, 33, 50, 18, 4, 49, 43, 28, 29, 47, 34, 32, 31, 25, 22, 20, 2, 35, 48, 24, 13, 10, 5, 45, 23, 16, 30, 14, 1, 19, 41, 26, 37, 17, 11} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Plot the solution :[font = text; inactive; preserveAspect; animationSpeed = 0; ] Having done all the CPU intensive work with the imported fortran routine, we can now use Mathematica's graphics capabilities to present the solution. ;[s] 3:0,0;89,1;100,0;150,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] route = positions[[Append[order,First[order]]]]; cities= {PointSize[0.03],Point/@positions}; Show[ Graphics[{cities,Line[route]}], AspectRatio->1 ]; :[font = postscript; PostScript; formatAsPostScript; output; inactive; preserveAspect; pictureLeft = 98; pictureWidth = 282; pictureHeight = 282; endGroup; endGroup; animationSpeed = 0; ] %! %%Creator: Mathematica %%AspectRatio: 1 MathPictureStart /Courier findfont 10 scalefont setfont % Scaling calculations 0.0238095 0.960061 0.0238095 0.953334 [ [ 0 0 0 0 ] [ 1 1 0 0 ] ] MathScale % Start of Graphics 1 setlinecap 1 setlinejoin newpath %%Object: Graphics [ ] 0 setdash 0 setgray gsave grestore 0 0 moveto 1 0 lineto 1 1 lineto 0 1 lineto closepath clip newpath gsave gsave gsave 0.03 setlinewidth 0.52976 0.87418 Mdot 0.38959 0.65873 Mdot 0.46544 0.449 Mdot 0.25902 0.36796 Mdot 0.21582 0.90183 Mdot 0.45104 0.49571 Mdot 0.68337 0.4509 Mdot 0.83026 0.20685 Mdot 0.46736 0.55291 Mdot 0.1515 0.97619 Mdot 0.97619 0.83414 Mdot 0.52592 0.39084 Mdot 0.09197 0.93138 Mdot 0.4376 0.92662 Mdot 0.58833 0.54814 Mdot 0.35599 0.78648 Mdot 0.96275 0.85512 Mdot 0.26478 0.4223 Mdot 0.54032 0.95617 Mdot 0.34831 0.62536 Mdot 0.63921 0.1058 Mdot 0.26575 0.64348 Mdot 0.33967 0.82938 Mdot 0.02573 0.90564 Mdot 0.1707 0.54719 Mdot 0.76498 0.73214 Mdot 0.70545 0.02762 Mdot 0.28879 0.14202 Mdot 0.22254 0.03144 Mdot 0.39055 0.86751 Mdot 0.11982 0.57674 Mdot 0.05741 0.51096 Mdot 0.4232 0.40324 Mdot 0.07085 0.45281 Mdot 0.21582 0.74739 Mdot 0.90227 0.60153 Mdot 0.78898 0.77504 Mdot 0.84754 0.56244 Mdot 0.75154 0.1201 Mdot 0.91571 0.33555 Mdot 0.60849 0.93234 Mdot 0.85618 0.59676 Mdot 0.37135 0.30123 Mdot 0.92915 0.20208 Mdot 0.30607 0.87609 Mdot 0.58257 0.18397 Mdot 0.10925 0.11056 Mdot 0.14478 0.81603 Mdot 0.35311 0.34413 Mdot 0.32143 0.4347 Mdot grestore grestore 0.004 setlinewidth 0.90227 0.60153 moveto 0.85618 0.59676 lineto 0.84754 0.56244 lineto 0.91571 0.33555 lineto 0.92915 0.20208 lineto 0.83026 0.20685 lineto 0.75154 0.1201 lineto 0.70545 0.02762 lineto 0.63921 0.1058 lineto 0.58257 0.18397 lineto 0.52592 0.39084 lineto 0.68337 0.4509 lineto 0.58833 0.54814 lineto 0.46736 0.55291 lineto 0.45104 0.49571 lineto 0.46544 0.449 lineto 0.4232 0.40324 lineto 0.32143 0.4347 lineto 0.26478 0.4223 lineto 0.25902 0.36796 lineto 0.35311 0.34413 lineto 0.37135 0.30123 lineto 0.28879 0.14202 lineto 0.22254 0.03144 lineto 0.10925 0.11056 lineto 0.07085 0.45281 lineto 0.05741 0.51096 lineto 0.11982 0.57674 lineto 0.1707 0.54719 lineto 0.26575 0.64348 lineto 0.34831 0.62536 lineto 0.38959 0.65873 lineto 0.21582 0.74739 lineto 0.14478 0.81603 lineto 0.02573 0.90564 lineto 0.09197 0.93138 lineto 0.1515 0.97619 lineto 0.21582 0.90183 lineto 0.30607 0.87609 lineto 0.33967 0.82938 lineto 0.35599 0.78648 lineto 0.39055 0.86751 lineto 0.4376 0.92662 lineto 0.52976 0.87418 lineto 0.54032 0.95617 lineto 0.60849 0.93234 lineto 0.76498 0.73214 lineto 0.78898 0.77504 lineto 0.96275 0.85512 lineto 0.97619 0.83414 lineto Mistroke 0.90227 0.60153 lineto Mfstroke grestore % End of Graphics MathPictureEnd :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Compute the path length :[font = text; inactive; preserveAspect; animationSpeed = 0; ] Note, due to round-off in the annealing algorithm, some computers may give a slightly different route than the one found above. Because the routine is based on random numbers it does not always find the optimal solution, but always something very close to it. :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Plus @@ Map[Sqrt[#.#]&, route-RotateLeft[route]] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; animationSpeed = 0; ] 5.740677969172233 ;[o] 5.74068 :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Solve the same problem but with a river :[font = text; inactive; preserveAspect; animationSpeed = 0; ] In this example we shall place a river along the line y=0.3 and impose a high penalty for crossing it. To do this we modify the cost function ALEN to include a penalty of lambda if the river is crossed, where the value of lambda can be adjusted to give the relative importance of path-length versus river-crossing. If lambda is high enough then the solution will cross the river a minimum number of times, ie. just twice. :[font = text; inactive; preserveAspect; animationSpeed = 0; ] It is possible to write the penalty function as a normal Mathematica function. Even though the external anneal routine is a fortran program it will, with InterCall, still be possible to call the MathematicaÐwritten penalty function. This ability to seamlessly link fortran or C code with Mathematica code is one of the features that makes InterCall so useful. ;[s] 7:0,0;57,1;68,0;195,1;206,0;288,1;299,0;360,-1; 2:4,13,9,Times,0,12,0,0,0;3,13,9,Times,2,12,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] penalty = Function[{x1,x2,y1,y2},Module[{lambda=1,m1,m2}, m1=If[y1<0.3,.5,-.5]; m2=If[y2<0.3,.5,-.5]; Sqrt[(x2-x1)^2+(y2-y1)^2]+lambda*(m2-m1)^2]]; order = anneal[x,y,$ALEN->penalty] :[font = output; output; inactive; preserveAspect; endGroup; ] {43, 49, 4, 34, 32, 31, 25, 18, 50, 33, 3, 12, 7, 15, 6, 9, 2, 20, 22, 35, 48, 24, 13, 10, 5, 45, 23, 16, 30, 14, 1, 19, 41, 26, 37, 17, 11, 36, 42, 38, 40, 44, 8, 39, 27, 21, 46, 29, 47, 28} ;[o] {43, 49, 4, 34, 32, 31, 25, 18, 50, 33, 3, 12, 7, 15, 6, 9, 2, 20, 22, 35, 48, 24, 13, 10, 5, 45, 23, 16, 30, 14, 1, 19, 41, 26, 37, 17, 11, 36, 42, 38, 40, 44, 8, 39, 27, 21, 46, 29, 47, 28} :[font = text; inactive; preserveAspect; endGroup; animationSpeed = 0; ] The above just looks like a normal Mathematica interaction, but of course underneath fortran code is being called (possibly on a remote computer) and also the penalty function is being remotely interpreted outside Mathematica. InterCall hides all these inner details from the user, who need only worry about writing Mathematica code. All other details are handled by InterCall. ;[s] 7:0,0;35,1;46,0;214,1;225,0;316,1;327,0;378,-1; 2:4,13,9,Times,0,12,0,0,0;3,13,9,Times,2,12,0,0,0; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Plot the solution with a river :[font = text; inactive; preserveAspect; animationSpeed = 0; ] As before we can plot the solution using Mathematica graphics: ;[s] 3:0,0;41,1;52,0;63,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] route = positions[[Append[order,First[order]]]]; river = {Thickness[0.03],GrayLevel[0.7],Line[{{0,0.3},{1,.3}}]}; cities= {PointSize[0.03],Point/@positions}; Show[ Graphics[{river,cities,Line[route]}], AspectRatio->1 ]; :[font = postscript; PostScript; formatAsPostScript; output; inactive; preserveAspect; pictureLeft = 98; pictureWidth = 282; pictureHeight = 282; endGroup; animationSpeed = 0; ] %! %%Creator: Mathematica %%AspectRatio: 1 MathPictureStart /Courier findfont 10 scalefont setfont % Scaling calculations 0.0238095 0.952381 0.0238095 0.953334 [ [ 0 0 0 0 ] [ 1 1 0 0 ] ] MathScale % Start of Graphics 1 setlinecap 1 setlinejoin newpath %%Object: Graphics [ ] 0 setdash 0 setgray gsave grestore 0 0 moveto 1 0 lineto 1 1 lineto 0 1 lineto closepath clip newpath gsave gsave 0.7 setgray 0.03 setlinewidth 0.02381 0.30981 moveto 0.97619 0.30981 lineto stroke grestore gsave gsave 0.03 setlinewidth 0.52571 0.87418 Mdot 0.38667 0.65873 Mdot 0.4619 0.449 Mdot 0.25714 0.36796 Mdot 0.21429 0.90183 Mdot 0.44762 0.49571 Mdot 0.6781 0.4509 Mdot 0.82381 0.20685 Mdot 0.46381 0.55291 Mdot 0.15048 0.97619 Mdot 0.96857 0.83414 Mdot 0.5219 0.39084 Mdot 0.09143 0.93138 Mdot 0.43429 0.92662 Mdot 0.58381 0.54814 Mdot 0.35333 0.78648 Mdot 0.95524 0.85512 Mdot 0.26286 0.4223 Mdot 0.53619 0.95617 Mdot 0.34571 0.62536 Mdot 0.63429 0.1058 Mdot 0.26381 0.64348 Mdot 0.33714 0.82938 Mdot 0.02571 0.90564 Mdot 0.16952 0.54719 Mdot 0.75905 0.73214 Mdot 0.7 0.02762 Mdot 0.28667 0.14202 Mdot 0.22095 0.03144 Mdot 0.38762 0.86751 Mdot 0.11905 0.57674 Mdot 0.05714 0.51096 Mdot 0.42 0.40324 Mdot 0.07048 0.45281 Mdot 0.21429 0.74739 Mdot 0.89524 0.60153 Mdot 0.78286 0.77504 Mdot 0.84095 0.56244 Mdot 0.74571 0.1201 Mdot 0.90857 0.33555 Mdot 0.60381 0.93234 Mdot 0.84952 0.59676 Mdot 0.36857 0.30123 Mdot 0.9219 0.20208 Mdot 0.30381 0.87609 Mdot 0.5781 0.18397 Mdot 0.10857 0.11056 Mdot 0.14381 0.81603 Mdot 0.35048 0.34413 Mdot 0.31905 0.4347 Mdot grestore grestore 0.004 setlinewidth 0.5781 0.18397 moveto 0.63429 0.1058 lineto 0.7 0.02762 lineto 0.74571 0.1201 lineto 0.82381 0.20685 lineto 0.9219 0.20208 lineto 0.90857 0.33555 lineto 0.84095 0.56244 lineto 0.84952 0.59676 lineto 0.89524 0.60153 lineto 0.96857 0.83414 lineto 0.95524 0.85512 lineto 0.78286 0.77504 lineto 0.75905 0.73214 lineto 0.60381 0.93234 lineto 0.53619 0.95617 lineto 0.52571 0.87418 lineto 0.43429 0.92662 lineto 0.38762 0.86751 lineto 0.35333 0.78648 lineto 0.33714 0.82938 lineto 0.30381 0.87609 lineto 0.21429 0.90183 lineto 0.15048 0.97619 lineto 0.09143 0.93138 lineto 0.02571 0.90564 lineto 0.14381 0.81603 lineto 0.21429 0.74739 lineto 0.38667 0.65873 lineto 0.34571 0.62536 lineto 0.26381 0.64348 lineto 0.16952 0.54719 lineto 0.11905 0.57674 lineto 0.05714 0.51096 lineto 0.07048 0.45281 lineto 0.25714 0.36796 lineto 0.26286 0.4223 lineto 0.31905 0.4347 lineto 0.44762 0.49571 lineto 0.46381 0.55291 lineto 0.58381 0.54814 lineto 0.6781 0.4509 lineto 0.5219 0.39084 lineto 0.4619 0.449 lineto 0.42 0.40324 lineto 0.35048 0.34413 lineto 0.36857 0.30123 lineto 0.28667 0.14202 lineto 0.10857 0.11056 lineto 0.22095 0.03144 lineto Mistroke 0.5781 0.18397 lineto Mfstroke grestore % End of Graphics MathPictureEnd :[font = text; inactive; preserveAspect; endGroup; animationSpeed = 0; ] To illustrate the difference in the solution with the river constraint added, you could select the two solutions and animate them. This will highlight the minor difference that is needed to avoid the river. :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Compute the path length of the above solution :[font = input; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Plus @@ Map[Sqrt[#.#]&, route-RotateLeft[route]] :[font = output; output; inactive; preserveAspect; endGroup; animationSpeed = 0; ] 5.901934589957522 ;[o] 5.90193 :[font = text; inactive; preserveAspect; endGroup; endGroup; animationSpeed = 0; ] The path length is slightly larger that before, as this is the price to pay if you don't like crossing rivers. :[font = section; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] Exercises :[font = subsubsection; inactive; preserveAspect; animationSpeed = 0; ] Modify the anneal.f program so that it can solve 3-dimensional problems. Use it to solve a "fruit picking problem". Plot (and maybe animate) your answer. :[font = subsubsection; inactive; preserveAspect; animationSpeed = 0; ] Try an animation problem by moving one of the cities around and seeing how the solution changes. Tip Ð because each iteration in the animation will be close to the previous solution, use the option $IORDER->order and also reduce the temperature $T to make the computation faster. This will require some experimentation. ;[s] 5:0,0;198,1;212,0;245,1;247,0;320,-1; 2:3,13,9,Times,1,12,0,0,0;2,13,10,Courier,1,12,0,0,0; :[font = subsubsection; inactive; preserveAspect; endGroup; animationSpeed = 0; ] Solve the 50 city problem by eye. An easy way to do this is to select the graphic cell showing the 50 cities and command-click on each city in the order you want. Then use copy (command-c), and paste (command-v) into an input cell. This will give you the co-ordinates, in order chosen. :[font = section; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] References :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] ANNEAL :[font = text; inactive; preserveAspect; animationSpeed = 0; ] The ANNEAL routine is an adaptation of a routine in the Numerical Recipes book by William H. Press, Brian P. Flannery, Saul A. Teukolsky and William T. Vetterling, 1987, 1992 New York, Cambridge University Press. Diskettes containing all the Numerial Recipes software can ordered from: :[font = text; inactive; preserveAspect; endGroup; ] Numerical Recipes Software or Customer Services Dept., Cambridge University, Cambridge University Press, 110 Midland Avenue, Edinburgh Building, Shaftesbury Road, Port Chester, NY 10573, Cambridge CB2 2RU, USA. UK. :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] InterCall :[font = text; inactive; preserveAspect; animationSpeed = 0; ] InterCall was developed by Terry Robb. For more information on InterCall send email to intercall_forward@wri.com. An order form is available via MathSource Ð email the line send 0202-587-0011 to mathsource@wri.com for details. Information can also be posted or faxed to Analytica below. ;[s] 9:0,0;87,1;112,0;173,1;191,0;195,1;213,0;270,2;279,0;287,-1; 3:5,13,9,Times,0,12,0,0,0;3,13,9,Times,1,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = text; inactive; preserveAspect; endGroup; endGroup; ] Analytica Fax: +61 (9) 386 5666 P.O. Box 522 Nedlands 6009 Perth WA AUSTRALIA. ;[s] 2:0,1;28,0;99,-1; 2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = section; inactive; Cclosed; pageBreak; preserveAspect; startGroup; animationSpeed = 0; ] Appendix :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; animationSpeed = 0; ] The modifications to the original anneal.f program :[font = text; inactive; preserveAspect; ] The external routine used in this Notebook, is a slightly modified version of the anneal subroutine given in the Numerical Recipes book. The routine has been converted to double precision and some control parameters have been passed as arguments and not hard-wired as in the original. :[font = text; inactive; preserveAspect; endGroup; endGroup; ] Specifically we have used real*8 statements for all real variables, and the original ALEN inline function has been passed as an argument to ANNEAL, TRNCST and REVCST. Workspace (named IWRK) for use by TRNSPT has also be passed as an argument to ANNEAL. Finally the number of T steps, which in the original was hard-wired to 100, has been changed to NTSTEPS and passed as an argument too. ^*)