Beginnings of a program to run approximation algorithms on instances of the traveling salesman problem.

tsp.lisp

;; (load "tsp.lisp")

(ql:quickload 'cl-ppcre)

(defstruct svg
  (points nil :type list)
  (polys nil :type list)
  (lines nil :type list))

(defstruct point
  (x 0.0 :type (or  single-float integer))
  (y 0.0 :type (or  single-float integer)))

(defun scale-point (pt s)
  (declare (point pt))
  (declare (number s))
  (make-point :x (* (point-x pt) s) :y (* (point-y pt) s)))

(defstruct line
  (pt1 nil :type point)
  (pt2 nil :type point))

(defun scale-line (ln s)
  (declare (line ln))
  (declare (number s))
  (make-line :pt1 (scale-point (line-pt1 ln) s) :pt2 (scale-point (line-pt2 ln) s)))
		   
(defstruct polygon
  (points nil :type list))

(defun add-point-coords (img x y)
  (declare (single-float x))
  (declare (single-float y))
  (setf (svg-points img) (cons (make-point :x x :y y) (svg-points img))))

(defun add-point (img pt)
  (setf (svg-points img) (cons pt (svg-points img))))

(defun add-line-coords (img x1 y1 x2 y2)
  (setf (svg-lines img) (cons (make-line :pt1 (make-point :x x1 :y y1) :pt2 (make-point :x x2 :y y2)) (svg-lines img))))

(defun add-line (img ln)
  (setf (svg-lines img) (cons ln (svg-lines img))))

(defun add-line-points (img pt1 pt2)
  (setf (svg-lines img) (cons (make-line :pt1 pt1 :pt2 pt2) (svg-lines img))))

(defun add-polygon-coords (img &rest pts)
  (setf
   (svg-polys img)
   (cons
	(make-polygon :points (mapcar #'(lambda (pt) (make-point :x (car pt) :y (cadr pt))) pts)) (svg-polys img))))

(defun to_file (img &key (stream t))
  (format stream "<?xml version=\"1.0\" standalone=\"no\"?> <!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\"> <svg width=\"100%\" height=\"100%\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">")
  (dolist (poly (svg-polys img))
	(format stream "<polygon points=\"")
	(dolist (pt (polygon-points poly))
	  (format stream "~a,~a " (point-x pt) (point-y pt)))
	(format stream "\" style=\"fill:#cccccc; stroke:#000000;stroke-width:1\"/>"))

  (dolist (line (svg-lines img))
	(let* ((p1 (line-pt1 line))
		   (p2 (line-pt2 line)))
			   
	(format stream "<line x1=\"~a\" y1=\"~a\" x2=\"~a\" y2=\"~a\" style=\"stroke:rgb(99,99,99);stroke-width:2\"/>" (point-x p1) (point-y p1) (point-x p2) (point-y p2))))
			
  (dolist (pt (svg-points img))
	(format stream "<circle cx=\"~a\" cy=\"~a\" r=\"3\" stroke=\"blue\" stroke-width=\"2\" fill=\"green\"/>" (point-x pt) (point-y pt)))
  (format stream "</svg>"))

(defun euclidean-distance (cities city1 city2)
  (let* ((pt1 (elt cities city1))
		 (pt2 (elt cities city2))
		 (x1 (point-x pt1))
		 (y1 (point-y pt1))
		 (x2 (point-x pt2))
		 (y2 (point-y pt2))
		 (xd (- x1 x2))
		 (yd (- y1 y2)))
	(declare (single-float x1))
	(declare (single-float y1))
	(declare (single-float x2))
	(declare (single-float y2))
	(declare (single-float xd))
	(declare (single-float yd))
	(sqrt (+ (* xd xd) (* yd yd)))))

(defun att-distance (cities city1 city2)
  (let* ((pt1 (elt cities city1))
		 (pt2 (elt cities city2))
		 (x1 (point-x pt1))
		 (y1 (point-y pt1))
		 (x2 (point-x pt2))
		 (y2 (point-y pt2))
		 (xd (- x1 x2))
		 (yd (- y1 y2))
		 (rij (sqrt (+ (* xd xd) (* yd yd))))
		 (tij (round rij)))
	(declare (single-float x1))
	(declare (single-float y1))
	(declare (single-float x2))
	(declare (single-float y2))
	(declare (single-float xd))
	(declare (single-float yd))
	(if (< tij rij)
		(+ tij 1)
	  tij)))

(defun total-distance (cities &key (distance-function #'euclidean-distance))
  (let ((distance 0.0)
		 (maxlen (- (length cities) 1)))
	(declare (single-float distance))
	(declare (integer maxlen))
	(loop for idx from 1 upto maxlen do
		  (setf distance (+ distance (funcall distance-function cities idx (- idx 1)))))
	(+ distance (funcall distance-function cities maxlen 0))))

(defun read-cities (fname)
  (let ((cities (make-array 1 :fill-pointer 0 :element-type 'point :adjustable t)))
	(with-open-file
	 (stream fname)
	 (loop for line = (read-line stream nil)
		   while line do
		   (let* ((string-coords (ppcre:split "\\s+" line))
				  (coords (mapcar #'read-from-string string-coords))
				  (pt (make-point :x (car coords) :y (cadr coords))))
			 (vector-push-extend pt cities))))
	cities))

(defun write-cities (cities fname)
  (with-open-file (stream fname :direction :output :if-exists :overwrite :if-does-not-exist :create)
				  (loop for city across cities do
						(format stream "~,8f ~,8f~%" (point-x city) (point-y city)))))

(defun cities-to-svg (cities &key (scale 1.0))
  (declare (single-float scale))
  (let ((img (make-svg))
		(maxlen (- (length cities) 1)))
		
	(declare (integer maxlen))
	(loop for city across cities do
		  (add-point img (scale-point city scale)))

	(loop for idx from 1 upto maxlen do
		  (add-line img (scale-line
						 (make-line
						  :pt1 (elt cities (- idx 1)) ; point at idx
						  :pt2 (elt cities idx)) scale))) ; point at idx + 1

	(add-line img (scale-line (make-line
						  :pt1 (elt cities maxlen)
						  :pt2 (elt cities 0)) scale))
	img))

(defun cities-to-svg-file (cities fname &key (scale 1.0))
  (with-open-file (stream fname :direction :output :if-exists :supersede :if-does-not-exist :create)
				  (to_file (cities-to-svg cities :scale scale) :stream stream)))

(defun swap-cities (cities i j)
  (if (not (equal i j))
	  (rotatef (aref cities i) (aref cities j))))

(defun three-dist (cities i &key (distance-function #'euclidean-distance))
  (let ((i1 (- i 1))
		(i2 i)
		(i3 (+ i 1))
		(max-len (- (length cities) 1)))
	(cond ((< i1 0 )
		   (+ (funcall distance-function cities max-len i2)
			  (funcall distance-function cities i2 i3)))
		  ((> i3 max-len )
		   (+ (funcall distance-function cities i1 i2)
			  (funcall distance-function cities i2 0)))
		  (t (+ (funcall distance-function cities i1 i2)
				(funcall distance-function cities i2 i3))))))

(defun two-opt (cities i j)
  (let ((old-total (total-3dist cities (list i j))))
	(swap-cities cities i j)
	(let ((new-total (total-3dist cities (list i j))))
	  (if (> new-total old-total)
		  (swap-cities cities i j)))))

(defun total-3dist (cities indices)
  (reduce #'+ (mapcar #'(lambda (x) (three-dist cities x)) indices)))

(defun do-3opt-perm (cities swaps)
  (dolist (swap swaps)
	(swap-cities cities (car swap) (cdr swap))))

;; [ 1 2 3 4 5 6 7 8 9 10 11 12]
;; [ 1 2 3 4 5 6 10 8 9 7 11 12]
;; [ 7 2 3 4 5 6 1 8 9 10 11 12]
;; [ 7 2 3 4 5 6 10 8 9 1 11 12]
;; [ 10 2 3 4 5 6 1 8 9 7 11 12]
;; [ 10 2 3 4 5 6 7 8 9 1 11 12]

(defun three-opt (cities i j k)
  (let ((all-combos (make-array 6 :initial-contents
								(list
								 (list nil )
								 (list (cons j k) )
								 (list (cons i j))
								 (list (cons i j) (cons j k))
								 (list (cons i k) )
								 (list (cons i k) (cons j k)))))
		(distances (make-array 6 :initial-contents '(0.0 0.0 0.0 0.0 0.0 0.0)))
		(3list (list i j k)))

	(dotimes (idx (length all-combos))
	  (do-3opt-perm cities (elt all-combos idx))
	  (setf (aref distances idx) (total-3dist cities 3list))
	  (do-3opt-perm cities (elt all-combos idx)))
	(let ((perm-to-use (position (reduce #'min distances) distances)))
	  (do-3opt-perm cities (elt all-combos perm-to-use)))))
	
(defun reduce-with-swaps (cities &key (max-swaps 500))
  (dotimes (i max-swaps)
	(let* ((maxlen (length cities))
		   (i1 (random maxlen))
		   (i2 (random maxlen))
		   (old-dist (total-distance cities)))
	  (swap-cities cities i1 i2)
	  (if (> (total-distance cities) old-dist)
		  (swap-cities cities i1 i2)))))

(defun reduce-with-2opt (cities &key (visualize nil))
  (let ((maxlen (length cities)))
  (dotimes (i maxlen)
	(dotimes (j maxlen)
	  (if (not (= i j))
		  (two-opt cities i j))
	  (if (and (= 0 (mod i 1000)) visualize)
		  (funcall visualize cities i))
	  ))))

(defun reduce-with-3opt (cities &key (visualize nil))
  (let ((maxlen (length cities)))
	(dotimes (i maxlen)
	  (dotimes (j maxlen)
		(dotimes (k maxlen)
		  (if (not (or (= i j) (= j k)))
			  (three-opt cities i j k))
		  (if (and (= 0 (mod i 1000)) visualize)
			  (funcall visualize cities i))
		  )))))

(defun reduce-with-3opt-random (cities &key (max-swaps 500) (visualize nil))
  (let ((maxlen (length cities)))
	(dotimes (i max-swaps)
	  (three-opt cities (random maxlen) (random maxlen) (random maxlen))
	  (if (and (= 0 (mod i 1000)) visualize)
		  (funcall visualize cities i))
	  )))

(defun reduce-with-2opt-random (cities &key (max-swaps 500) (visualize nil))
  (let ((maxlen (length cities)))
  (dotimes (i max-swaps)
	(let ((ri (random maxlen))
		  (rj (random maxlen)))
	  (two-opt cities ri rj)
	  (if (and (= 0 (mod i 1000)) visualize)
		  (funcall visualize cities i))
	  ))))

(defun randomize-cities (cities &key (max-swaps 500) (visualize nil))
  (dotimes (i max-swaps)
	(swap-cities cities (random (length cities)) (random (length cities)))
	(if (and (= 0 (mod i 1000)) visualize)
		  (funcall visualize cities i))
	))

(defun show-cities (cities iteration)
  (cities-to-svg-file cities (format nil "images/cities~a.svg" iteration) :scale 0.1))