; the Travelling Salesman Problem (https://en.wikipedia.org/wiki/Travelling_salesman_problem)
; a Simulated Annealing example (combinatorial minimisation), by RTFC (22 Feb 2016)
; Adapted from Press et al. Numerical Recipes, 2nd ed., pp.438-443.

; Note that the algorithm converges on A local minimum (in terms of the
; user-defined cost-function(s)), which is not necessarily THE global minimum.
; Note also that the search path, duration, and final result may differ from run to run (city coordinates change every time).
; Several parameters and weights can be tweaked to adjust this.

#include <Array.au3>
#include <GUIConstants.au3>
#include <GDIplus.au3>
#include <Math.au3>

Global Const $GUIwidth=600
Global Const $GUIheight=400
Global Const $citysize=10
Global Const $halfcitysize=$citysize/2

$hwnd=GUICreate("Travelling Salesman Problem (using Simulated Annealing), by RTFC",$GUIwidth,$GUIheight)
GUISetState()

; Load GDI+ resources
_GDIPlus_Startup()
Global $graphics=_GDIPlus_GraphicsCreateFromHWND($hwnd)
Global $bitmap=_GDIPlus_BitmapCreateFromGraphics($GUIwidth,$GUIheight,$graphics)
Global $backbuffer=_GDIPlus_ImageGetGraphicsContext($bitmap)
Global $hPenLine = _GDIPlus_PenCreate(0xFF0000FF, 4)
Global $hPenCircle = _GDIPlus_PenCreate(0xFFFF0000, 4)


; Simulated Annealing vars
Global $temperat,$path,$kk,$nswap,$nswapstep,$cost,$tempstep,$absimp,$lowestcost

; initialise
SRandom(@SEC+@AutoItPID); initialise randomising seed
$verbose=True			; T: write regular progress updates to console

; define cities
Global $ncity=20
if $ncity<10 then Exit		; we need something to work with

Global $maxcity=$ncity*50
Global $n[7],$xx[7],$yy[7]	; base-1 indexed
Global $x[$ncity+1],$y[$ncity+1],$iorder[$ncity+1],$jorder[$maxcity+1]	; base-1 indexed

;______START OF ANNEALING ROUTINE____________

$nover=100*$ncity		; maximum number of paths at any temperature
$nlimit=10*$ncity		; maximum number of successful path changes before continuing
$nwrite=Int($nover/5)	; default status update interval if verbose=.t.
$tempsteps=100			; number of temperature steps to try
$tfactor=0.90			; annealing schedule: temperature is reduced by this factor after each step
$path=0

While True
	$temperat=0.5		; initial temperature; smaller = more aggressive + more myopic search
	$absimp=0			; counter
	$nswapstepzero=0	; counter

	; prep the buffers
	For $cc=1 to 6
		$n[$cc]=0
	Next

	For $cc=1 to $ncity
		$x[$cc]=Random()
		$y[$cc]=Random()
		$iorder[$cc]=$cc
		$jorder[$cc]=0
	Next

	For $cc=$ncity+1 to $maxcity
		$jorder[$cc]=0
	Next

	; prep the cost vars
	_CalcPathLength()
	$initcost=$path
	$lowestcost=$path

	; main loop starts here
	For $tempstep=1 to $tempsteps		; try up to N temperature steps
		$nswap=0
		$nswapstep=0

		For $kk=1 to $nover
			$n[1]=Random(1,$ncity,1)   ; choose beginning and end of segment
			$n[2]=Random(1,$ncity,1)
			If $n[2]>=$n[1] Then $n[2]=1+Mod(($n[2]+1),$ncity)

			; count number of cities not on segment ($nn)
			$nn=1+Mod(($n[1]-$n[2]+$ncity-1),$ncity)
			If $nn<3 then ContinueLoop

			; decide whether to do a segment transport or reversal (equal chances)
			$doTransport=(Random()<=0.5)
			If $doTransport Then		; try a segment transport
				$n[3]=$n[2]+Int(Abs($nn-2)*Random())+1
				$n[3]=1+Mod($n[3]-1,$ncity)	; transport to a location not on the path
				_TransportCost()

			Else		; try a reversal instead
				_ReversalCost()
			EndIf

			; Listen to the wind, talk to the trees...
			Switch _AskOracle()
				Case True
					$nswap+=1
					$path+=$cost
					if $doTransport Then
						_DoTransport()
					Else
						_DoReversal()
					EndIf

					If $lowestcost>$path Then
						$nswapstep+=1
						$absimp+=1
						$lowestcost=$path
					EndIf

					If $nswap>=$nlimit then ExitLoop
			EndSwitch
		Next

		If $verbose Then _ScreenOut()
		If $nswapstep=0 then $nswapstepzero+=1
		If $nswapstepzero=30 then ExitLoop		; no more improvements in the last N temperature steps

		; reduce temperature = likelihood of following a trajectory away from the nearest LOCAL optimum (in the hope of getting nearer to the GLOBAL optimum)
		$temperat*=$tfactor
	Next

	; show final result
	if Msgbox($MB_OKCANCEL,"Simulated Annealing Test Result","Shortest Path Length: " & $lowestcost)=$IDCANCEL then Exit

WEnd

_GDIclose()
Exit


Func _AskOracle()

	If $cost<0 Then
		Return True
	Else		; this is where all the magic happens!
		Return (random()<Exp(-($cost/$temperat)))
	Endif

EndFunc


Func _CalcPathLength()

  $path=0	; Global var
  Local $cc,$index1,$index2
	For $cc=1 To $ncity-1
		$index1=$iorder[$cc]
		$index2=$iorder[$cc+1]
		$path+=_Distance($x[$index1],$x[$index2],$y[$index1],$y[$index2])
	Next
	$index1=$iorder[$ncity]      ; close the loop by tying path ends together
	$index2=$iorder[1]
	$path+=_Distance($x[$index1],$x[$index2],$y[$index1],$y[$index2])

EndFunc


Func _Distance($x1,$x2,$y1,$y2)

	Return Sqrt((($x1-$x2)^2)+(($y1-$y2)^2))
EndFunc


Func _ReversalCost()

	$n[3]=1+Mod(($n[1]+$ncity-2),$ncity)	; find the city before n[1]...
	$n[4]=1+Mod($n[2],$ncity)				; and after n[2]

	; find coordinates of the four cities
	local $ii,$jj,$dex
	For $jj=1 to 4
		$dex=$n[$jj]
		$ii=$iorder[$dex]
		$xx[$jj]=$x[$ii]
		$yy[$jj]=$y[$ii]
	Next

	; calculate the cost of disconnecting the segment at both ends and reconnecting in the opposite order
	$cost = -_Distance($xx[1],$xx[3],$yy[1],$yy[3])
	$cost-=_Distance($xx[2],$xx[4],$yy[2],$yy[4])
	$cost+=_Distance($xx[1],$xx[4],$yy[1],$yy[4])
	$cost+=_Distance($xx[2],$xx[3],$yy[2],$yy[3])

EndFunc


Func _DoReversal()

	Local $mm=(1+Mod(($n[2]-$n[1]+$ncity),$ncity))*.5
	local $ii,$jj,$ll,$itmp
	For $jj=1 to $mm
		$ii=1+Mod(($n[1]+$jj-2),$ncity)
		$ll=1+Mod(($n[2]-$jj+$ncity),$ncity)

		$itmp=$iorder[$ii]
		$iorder[$ii]=$iorder[$ll]
		$iorder[$ll]=$itmp
	Next

EndFunc


Func _TransportCost()

	$n[4]=1+Mod($n[3],$ncity)				; find the city following n[3]
	$n[5]=1+Mod(($n[1]+$ncity-2),$ncity)	; find the city before n[1]...
	$n[6]=1+Mod($n[2],$ncity) 				; and after n[2]

	; find coordinates of the six cities
	local $jj,$dex
	For $jj=1 to 6
		$dex=$n[$jj]
		$ii=$iorder[$dex]
		$xx[$jj]=$x[$ii]
		$yy[$jj]=$y[$ii]
	Next

	; calculate the cost of disconnecting the path segment n1->n2, opening a space between n3-n4, connecting the segment in the space, and conencting n5 to n6
	$cost=-_Distance($xx[2],$xx[6],$yy[2],$yy[6])
	$cost-=_Distance($xx[1],$xx[5],$yy[1],$yy[5])
	$cost-=_Distance($xx[3],$xx[4],$yy[3],$yy[4])
	$cost+=_Distance($xx[1],$xx[3],$yy[1],$yy[3])
	$cost+=_Distance($xx[2],$xx[4],$yy[2],$yy[4])
	$cost+=_Distance($xx[5],$xx[6],$yy[5],$yy[6])

EndFunc


Func _DoTransport()

	Local $m1,$m2,$m3,$pp
	$m1=1+Mod(($n[2]-$n[1]+$ncity),$ncity)	; find # cities from n1->n2
	$m2=1+Mod(($n[5]-$n[4]+$ncity),$ncity)	; find # cities from n4->n5
	$m3=1+Mod(($n[3]-$n[6]+$ncity),$ncity)	; find # cities from n6->n3

	Local $mm=1
	For $pp=1 to $m1
		$jj=1+Mod($pp+$n[1]-2,$ncity)
		$jorder[$mm]=$iorder[$jj]
		$mm+=1
	Next

	For $pp=1 to $m2
		$jj=1+Mod($pp+$n[4]-2,$ncity)
		$jorder[$mm]=$iorder[$jj]
		$mm+=1
	Next

	For $pp=1 to $m3
		$jj=1+Mod($pp+$n[6]-2,$ncity)
		$jorder[$mm]=$iorder[$jj]
		$mm+=1
	Next

	for $pp=1 to $ncity
		$iorder[$pp]=$jorder[$pp]
	Next

EndFunc


Func _ScreenOut()

	ConsoleWrite("Simulated Annealing. Initial Path length: " & $initcost & @CRLF)
	ConsoleWrite("Step: " & $tempstep & " of " & $tempsteps & "; Temperature: " & $temperat & @CRLF)
	ConsoleWrite("Executed Swaps: " & $nswap & "; shortest Path length: " & $lowestcost & @CRLF)
	ConsoleWrite("Total Improvements: " & $absimp & "; Improvements this step: " & $nswapstep & @CRLF  & @CRLF)
	_DrawTSroute()

EndFunc


func _DrawTSroute()

	_GDIPlus_GraphicsClear($backbuffer)

	For $cc=1 To $ncity-1
		$index1=$iorder[$cc]
		$index2=$iorder[$cc+1]
		_DrawLine($x[$index1],$x[$index2],$y[$index1],$y[$index2])
	Next
	$index1=$iorder[$ncity]      ; close the loop by tying path ends together
	$index2=$iorder[1]
	_DrawLine($x[$index1],$x[$index2],$y[$index1],$y[$index2])
    _GDIPlus_GraphicsDrawImageRect($graphics,$bitmap,0,0,$GUIwidth,$GUIheight)

EndFunc


Func _DrawLine($x1,$x2,$y1,$y2)

	local $x1scaled=Round($GUIwidth*$x1,0)
	Local $y1scaled=Round($GUIheight*$y1,0)
	Local $x2scaled=Round($GUIwidth*$x2,0)
	Local $y2scaled=Round($GUIheight*$y2,0)

	_GDIPlus_GraphicsDrawArc($backbuffer, $x1scaled-$halfcitysize, $y1scaled-$halfcitysize, $citysize, $citysize, 0, 360, $hPenCircle)
    _GDIPlus_GraphicsDrawLine($backbuffer, $x1scaled, $y1scaled, $x2scaled, $y2scaled, $hPenLine)

EndFunc


Func _GDIclose()

	_GDIPlus_PenDispose($hPenLine)
	_GDIPlus_PenDispose($hPenCircle)
    _GDIPlus_GraphicsDispose($backbuffer)
    _GDIPlus_BitmapDispose($bitmap)
    _GDIPlus_GraphicsDispose($graphics)
    _GDIPlus_Shutdown()

EndFunc