share/graphs/wiener_index.lisp

   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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, 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")))))