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