1;;;
2;;;  GRAPHS - graph theory package for Maxima
3;;;
4;;;  Copyright (C) 2008 Andrej Vodopivec <andrej.vodopivec@gmail.com>
5;;;
6;;;  This program is free software; you can redistribute it and/or modify
7;;;  it under the terms of the GNU General Public License as published by
8;;;  the Free Software Foundation; either version 2 of the License, or
9;;;  (at your option) any later version.
10;;;
11;;;  This program is distributed in the hope that it will be useful,
12;;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
13;;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14;;;  GNU General Public License for more details.
15;;;
16;;;  You should have received a copy of the GNU General Public License
17;;;  along with this program; if not, write to the Free Software
18;;;  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19;;;
20
21
22;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
23;;;
24;;; This file contains programs for computing the
25;;; Wiener index of a graph. It includes the algorithms for
26;;; ordinary and weighted Wiener index computation.
27;;;
28;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
29
30;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
31;;;
32;;; FLOYD-WARSHALL algorithm for all-pairs shortest paths
33;;;
34;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
35
36(in-package :maxima)
37
38(defun find-option (opt options &optional default)
39  (dolist (o options)
40    (cond ((eq o opt) (return-from find-option t))
41	  ((and (listp o) (eq (cadr o) opt))
42	   (return-from find-option (caddr o)))))
43  default)
44
45(defmfun $floyd_warshall (g &rest options)
46  (require-graph-or-digraph 1 'floyd_warshall g)
47  (let* ((n (graph-order g))
48	 (m (make-array (list n n)))
49	 (my-inf 1)
50	 (options (cons '((mlist simp)) options))
51	 (weighted (find-option '$weighted options t))
52	 (vertices (cdr (find-option '$vertices options)))
53	 (mat ($zeromatrix n n)))
54
55    (unless vertices
56      (setq vertices (vertices g)))
57
58    (if weighted
59	(dolist (e (cdr ($edges g)))
60	  (setq my-inf (m+ my-inf ($abs ($get_edge_weight e g 1)))))
61	(setq my-inf (1+ (graph-order g))))
62
63    ;; setup the array
64    (dotimes (i n)
65      (dotimes (j n)
66	(if (/= i j)
67	    (setf (aref m i j)
68		  (if weighted
69		      ($get_edge_weight `((mlist simp) ,(nth i vertices) ,(nth j vertices)) g
70					1 my-inf)
71		      (if (member (nth i vertices) (neighbors (nth j vertices) g)) 1 my-inf)))
72	    (setf (aref m i j) 0))))
73
74    ;; compute the distances
75    (dotimes (k n)
76      (dotimes (i n)
77	(dotimes (j n)
78	  (when (eq (mlsp (m+ (aref m i k) (aref m k j)) (aref m i j)) t)
79	    (setf (aref m i j) (m+ (aref m i k) (aref m k j)))))))
80
81    ;; check for negative cycles
82    (dotimes (k n)
83      (when (eq (mlsp (aref m k k) 0) t)
84	($error "Graph contains a negative cycle.")))
85
86    ;; fill the matrix
87    (dotimes (i n)
88      (dotimes (j n)
89	(setf (nth (1+ j) (nth (1+ i) mat))
90	      (if (equal (aref m i j) my-inf) '$inf (aref m i j)))))
91
92    mat))
93
94;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
95;;;
96;;; JOHNSON's algorithm for all pairs shortest path
97;;;
98;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
99
100(defun bellman-ford (s g)
101  (let ((d (make-hash-table))
102	(edges (if (graph-p g) (append (edges g) (mapcar #'reverse (edges g))) (edges g)))
103	(my-inf 1)
104	(prev (make-hash-table)))
105
106    (dolist (e (cdr ($edges g)))
107      (setf my-inf (m+ my-inf ($abs ($get_edge_weight e g 1)))))
108
109    ;; initialize distances
110    (dolist (v (vertices g))
111      (if (= s v)
112	  (setf (gethash v d) 0)
113	  (setf (gethash v d) my-inf)))
114
115    ;; relax edges
116    (dotimes (i (1- (length (vertices g))))
117      (dolist (e edges)
118	(let* ((u (first e))
119	       (v (second e))
120	       (nd (m+ (gethash u d) ($get_edge_weight `((mlist simp) ,@e) g))))
121	  (when (and (not (equal (gethash u d) my-inf))
122		     (eq (mgrp (gethash v d) nd) t))
123	    (setf (gethash v d) nd)
124	    (setf (gethash v prev) u)))))
125
126    ;; check for negative cycles
127    (dolist (e edges)
128      (let ((u (first e))
129	    (v (second e)))
130	(when (eq (mgrp (gethash v d)
131			(m+ (gethash u d) ($get_edge_weight `((mlist simp) ,@e) g)))
132		  t)
133	  ($error "Graph contains a negative cycle."))))
134
135    (values d prev)))
136
137(defmfun $johnson (g &rest options)
138  (let* ((h ($copy_graph g))
139	 (n (graph-order g))
140	 (options (cons '((mlist simp)) options))
141	 (weighted (find-option '$weighted options t))
142	 (vertices (cdr (find-option '$vertices options)))
143	 (m ($zeromatrix n n))
144	 nv)
145
146    (unless vertices
147      (setq vertices (vertices g)))
148    (setq nv (1+ (apply #'max vertices)))
149
150    (when weighted
151      (dolist (e (cdr ($edges g)))
152	($set_edge_weight e ($get_edge_weight e g) h)))
153
154    ;; add a new vertex
155    ($add_vertex nv h)
156    (dolist (v vertices)
157      ($add_edge `((mlist simp) ,nv ,v) h)
158      ($set_edge_weight `((mlist simp) ,nv ,v) 0 h))
159
160    ;; run the bellman-ford algorithm
161    (multiple-value-bind (d prev)
162	(bellman-ford nv h)
163      (declare (ignore prev))
164
165      ;; re-weight the edges
166      (dolist (e (cdr ($edges g)))
167	(let ((nw (m+ (if weighted ($get_edge_weight e g) 1)
168		      (gethash ($first e) d)
169		      (m- (gethash ($second e) d)))))
170	($set_edge_weight e nw h)))
171      ($remove_vertex nv h)
172
173      ;; run the dijkstra's algorithm for each vertex
174      (dotimes (i n)
175	(multiple-value-bind (dd pd)
176	    (dijkstra (nth i vertices) nv h)
177	  (declare (ignore pd))
178	  (dotimes (j n)
179	    (when (/= i j)
180	      (setf (nth (1+ j) (nth (1+ i) m))
181		    (if (eq '$inf (gethash (nth j vertices) dd))
182			'$inf
183			(m+ (gethash (nth j vertices) dd)
184			    (m- (gethash (nth i vertices) d))
185			    (gethash (nth j vertices) d))))))))
186
187      ;; return the matrix m
188      m)))
189
190;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
191;;;
192;;; JUVAN-MOHAR algorithm
193;;;
194;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
195
196(defun wiener-index (g)
197  (let ((n (length (vertices g)))
198	(q (make-hash-table))
199	(j 0) (k 1) (w 0)
200	(visited (make-hash-table))
201	(d (make-hash-table :test #'equal))
202	(s (make-hash-table))
203	(c (make-hash-table :test #'equal))
204	(we (make-hash-table :test #'equal))
205	(wi 0))
206
207    ;; initialize array we
208    (dolist (e (edges g))
209      (setf (gethash e we) 0))
210
211    ;; initialize arrays c and d
212    (dolist (v (vertices g))
213      (dolist (u (vertices g))
214	(setf (gethash (list v u) c) 0)
215	(setf (gethash (list v u) d) 0)))
216
217    ;; For each vertex i
218    (dolist (v (vertices g))
219
220      ;; initialize c
221      (setf (gethash (list v v) c) 1)
222
223      ;; we do a bfs search
224      (dolist (i (vertices g))
225	(setf (gethash i visited) nil))
226      (setf (gethash 1 q) v
227	    (gethash v visited) t)
228      (setq j 0 k 1)
229      (loop while (< j k) do
230	   (incf j)
231	   (setq w (gethash j q))
232	   (dolist (u (neighbors w g))
233	     (unless (gethash u visited)
234	       (setf (gethash (list v u) d) (1+ (gethash (list v w) d)))
235	       (setf (gethash u visited) t)
236	       (incf k)
237	       (setf (gethash k q) u))
238	     (when (> (gethash (list v u) d) (gethash (list v w) d))
239	       (incf (gethash (list v u) c) (gethash (list v w) c)))))
240
241      ;; and visit vertices in reverse order
242      (loop for j from n downto 1 do
243	   (let ((w (gethash j q)))
244	     (setf (gethash w s) 0)
245	     (dolist (u (neighbors w g))
246	       (when (< (gethash (list v w) d)
247			(gethash (list v u) d))
248		 (let ((x (* (/ (gethash (list v w) c)
249				(gethash (list v u) c))
250			     (1+ (gethash u s)))))
251		   (incf (gethash (list (min u w) (max u w)) we) (/ x 2))
252		   (incf (gethash w s) x)))))) )
253
254    ;; Compute the Wiener index
255    (dolist (e (edges g))
256      (incf wi (gethash e we)))
257
258    wi))
259
260(defmfun $wiener_index (g &rest options)
261  (require-graph 1 'wiener_index g)
262  (unless ($is_connected g)
263    ($error "`wiener_index': input graph is not connected"))
264  (let* ((weighted (find-option '$weighted options nil))
265	 (algorithm (find-option '$algorithm options (if weighted '$floyd_warshall '$juvan_mohar))))
266    (case algorithm
267      ($juvan_mohar (wiener-index g))
268      ($johnson
269       (m// ($xreduce "+" ($flatten ($args ($johnson g `((mlist simp) $weighted ,weighted))))) 2))
270      ($floyd_warshall
271       (m// ($xreduce "+" ($flatten ($args ($floyd_warshall g `((mlist simp) $weighted ,weighted))))) 2))
272      (t ($error "Unknown algorithm for WIENER_INDEX")))))
273