share/graphs/matching.lisp

   1 ;;;
   2 ;;;  GRAPHS - graph theory package for Maxima
   3 ;;;
   4 ;;;  Copyright (C) 2007 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 (in-package :maxima)
  22
  23 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  24 ;;;
  25 ;;; matching algorithms
  26 ;;;
  27 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  28
  29
  30 (defmfun $max_matching (gr)
  31   (require-graph 'maximum_matching 1 gr)
  32   (let ((partition (cdr ($bipartition gr))))
  33     (if (null partition)
  34         (maximum-matching-general gr)
  35         (maximum-matching-bipartite gr (cdr (car partition)) (cdr (cadr partition))))))
  36
  37 ;;;;;;;;;;
  38 ;;
  39 ;; the bipartite case
  40 ;;
  41
  42 (defun maximum-matching-bipartite (gr a b &optional cover)
  43   (let ((matching (make-hash-table))
  44         (done))
  45
  46     ;; greedy matching
  47     (loop for e in (edges gr) do
  48          (let ((u (first e))
  49                (v (second e)))
  50            (when (and (null (gethash u matching))
  51                       (null (gethash v matching)))
  52              (setf (gethash u matching) v)
  53              (setf (gethash v matching) u))))
  54
  55     ;; augment the matching
  56     (loop while (not done) do
  57          (setq done t)
  58          (let ((au) (bu))
  59
  60            ;; find unmatched vertices
  61            (loop for u in a do
  62                 (when (null (gethash u matching))
  63                   (push u au)))
  64            (loop for u in b do
  65                 (when (null (gethash u matching))
  66                   (push u bu)))
  67
  68            ;; find augmenting path
  69            (when (not (or (null au)
  70                           (null bu)))
  71              (let ((prev (make-hash-table))
  72                    (new1 au) (new2) (x))
  73                (loop while (not (null new1)) do
  74
  75                   ;; add edges in M^c
  76                     (setq new2 ())
  77                     (loop for u in new1 do
  78                          (loop for v in (neighbors u gr) do
  79                               (when (and (null (gethash v prev))
  80                                          (or (null (gethash v matching))
  81                                              (not (= (gethash v matching) u))))
  82                                 (setf (gethash v prev) u)
  83                                 (push v new2))))
  84                     (setq new1 ())
  85
  86                   ;; add edges in M
  87                     (loop for v in new2 do
  88                          (let ((u (gethash v matching)))
  89                            (unless (null u)
  90                              (push u new1)
  91                              (setf (gethash u prev) v)))))
  92
  93                ;; check for augmenting path
  94                (loop for v in bu while (null x) do
  95                     (unless (null (gethash v prev))
  96                       (setq x v)))
  97
  98                ;; extend the matching
  99                (unless (null x)
 100                  (setq done nil)
 101                  (loop while (not (null x)) do
 102                       (let ((u x)
 103                             (v (gethash x prev)))
 104                         (setf (gethash u matching) v)
 105                         (setf (gethash v matching) u)
 106                         (setq x (gethash v prev))))) )) ))
 107     (if (null cover)
 108         ;; we want the matching
 109         (let ((mmatching ()))
 110           (maphash #'(lambda (u v) (if (< u v)
 111                                        (setq mmatching (cons `((mlist simp) ,u ,v) mmatching))))
 112                    matching)
 113           (cons '(mlist simp) mmatching))
 114         ;; we want the vertex cover
 115         (get-cover-from-matching gr a matching)) ))
 116
 117 (defun get-cover-from-matching (gr a matching)
 118   (let ((au) (cov))
 119     (loop for x in (cdr a) do
 120          (when (null (gethash x matching))
 121            (push x au)))
 122     (if (null au)
 123         `((mlist simp) ,@(sort (cdr a) #'<))
 124
 125         ;; construct the Hungarian tree
 126         (let ((prev (make-hash-table))
 127               (new1 au) (new2))
 128           (loop while (not (null new1)) do
 129              ;; add edges in M^c
 130                (setq new2 ())
 131                (loop for u in new1 do
 132                     (loop for v in (neighbors u gr) do
 133                          (when (and (null (gethash v prev))
 134                                     (or (null (gethash v matching))
 135                                         (not (= (gethash v matching) u))))
 136                            (setf (gethash v prev) u)
 137                            (push v new2))))
 138                      (setq new1 ())
 139
 140              ;; add edges in M
 141                (loop for v in new2 do
 142                     (let ((u (gethash v matching)))
 143                       (unless (null u)
 144                         (push u new1)
 145                         (setf (gethash u prev) v)))))
 146           (maphash #'(lambda (u v)
 147                        (when (member v a)
 148                          (if (null (gethash v prev))
 149                              (push v cov)
 150                              (push u cov))))
 151                    matching)
 152           `((mlist simp) ,@(sort cov #'<))))))
 153
 154 ;;;;;;;;;;;
 155 ;;
 156 ;; set partition using hash-tales
 157 ;;
 158
 159 (defun main-vertex (v sp)
 160   (let ((m v))
 161     (loop while (not (= m (gethash m sp))) do
 162          (setq m (gethash m sp)))
 163     (setf (gethash v sp) m)
 164     m))
 165
 166 (defun join-sets (v u sp)
 167   (setf (gethash (main-vertex u sp) sp)
 168         (main-vertex v sp)))
 169
 170 (defun set-main-vertex (m v sp)
 171   (setf (gethash (main-vertex v sp) sp) m)
 172   (setf (gethash m sp) m))
 173
 174 ;;;;;;;;;;;
 175 ;;
 176 ;; the general case
 177 ;;
 178
 179 (defvar *matching-state*)
 180 (defvar *matching-matching*)
 181 (defvar *matching-bridge*)
 182 (defvar *matching-prev*)
 183 (defvar *matching-path*)
 184
 185 (defun maximum-matching-general (gr)
 186   (let ((*matching-matching* (make-hash-table))
 187         (done))
 188
 189     ;; greedy matching
 190     (loop for e in (edges gr) do
 191          (let ((u (first e))
 192                (v (second e)))
 193            (when (and (null (gethash u *matching-matching*))
 194                       (null (gethash v *matching-matching*)))
 195              (setf (gethash u *matching-matching*) v)
 196              (setf (gethash v *matching-matching*) u))))
 197
 198     ;; augment the matching
 199     (loop while (not done) do
 200          (setq done t)
 201          (let ((unmatched-vertices)
 202                (dfsnum)
 203                (dfsnum-curr)
 204                (vertex-partition)
 205                (*matching-state*)
 206                (*matching-bridge* (make-hash-table))
 207                (active-edges)
 208                (*matching-prev* (make-hash-table))
 209                (x))
 210
 211            ;; find unmatched vertices
 212            (loop for v in (vertices gr) do
 213                 (when (null (gethash v *matching-matching*))
 214                   (push v unmatched-vertices)))
 215
 216            (loop for first-vertex in unmatched-vertices while (null x) do
 217                 (setq dfsnum (make-hash-table))
 218                 (setq dfsnum-curr 0)
 219                 (setq *matching-state* (make-hash-table))
 220                 (setq *matching-bridge* (make-hash-table))
 221                 (setq vertex-partition (make-hash-table))
 222                 (setq *matching-state* (make-hash-table))
 223                 (setf (gethash first-vertex *matching-state*) 'A)
 224                 (setq active-edges ())
 225
 226                 (loop for v in (vertices gr) do
 227                      (setf (gethash v vertex-partition) v))
 228                 (loop for u in (neighbors first-vertex gr) do
 229                      (push (list first-vertex u) active-edges))
 230
 231
 232               ;; find augmenting path
 233                 (loop while (and (not (null active-edges)) (null x)) do
 234                      (let* ((e (pop active-edges))
 235                             (v (first e))
 236                             (w (second e))
 237                             (v-main (main-vertex v vertex-partition))
 238                             (w-main (main-vertex w vertex-partition)))
 239
 240                        (when (null (gethash v-main dfsnum))
 241                          (setf (gethash v-main dfsnum) dfsnum-curr)
 242                          (incf dfsnum-curr))
 243
 244                        (cond
 245
 246                          ;; we found an agugmented path
 247                          ((and (member w-main unmatched-vertices)
 248                                (not (= w-main first-vertex)))
 249                           (setf (gethash w-main *matching-prev*) v)
 250                           (setq x w-main))
 251
 252                          ;; we didn't visit w yet
 253                          ((null (gethash w-main *matching-state*))
 254                           (setf (gethash w-main *matching-state*) 'B)
 255                           (let ((z (gethash w-main *matching-matching*)))
 256                             (setf (gethash z *matching-state*) 'A)
 257                             (loop for x in (neighbors z gr) do
 258                                  (unless (= x w-main)
 259                                    (push (list z x) active-edges))))
 260                           (setf (gethash w-main *matching-prev*) v))
 261
 262                          ;; found a blossom
 263                          ((and (not (= w-main v-main)) ;; in the same blossom already
 264                                (eql (gethash w-main *matching-state*) 'A))
 265
 266                           (let ((fst) (lst) (b))
 267                             (if (> (gethash w-main dfsnum)
 268                                    (gethash v-main dfsnum))
 269                                 (setq fst w-main
 270                                       lst v-main
 271                                       b (list w v))
 272                                 (setq fst v-main
 273                                       lst w-main
 274                                       b (list v w)))
 275
 276                             (let ((tmp-vrt fst))
 277                               (loop while (not (= tmp-vrt lst)) do
 278                                    (setf (gethash tmp-vrt *matching-bridge*) b)
 279                                    (join-sets lst tmp-vrt vertex-partition)
 280                                    (if (eql (gethash tmp-vrt *matching-state*) 'A)
 281
 282                                        (progn
 283                                          (setq tmp-vrt (main-vertex (gethash tmp-vrt *matching-matching*)
 284                                                                     vertex-partition))
 285                                          (loop for u in (neighbors tmp-vrt gr) do
 286                                               (setq u (main-vertex u vertex-partition))
 287                                               (push (list tmp-vrt u) active-edges)))
 288
 289                                        (progn
 290                                          (setq tmp-vrt (main-vertex (gethash tmp-vrt *matching-prev*)
 291                                                                     vertex-partition))))))) )
 292                          )))
 293
 294               ;; augment the path
 295                 (unless (null x)
 296                   (setq done nil)
 297                   (let ((*matching-path* (list x)))
 298                     (find-augmenting-path first-vertex (gethash x *matching-prev*))
 299                     (loop while *matching-path* do
 300                          (let ((a (car *matching-path*))
 301                                (b (cadr *matching-path*)))
 302                            (setf (gethash a *matching-matching*) b)
 303                            (setf (gethash b *matching-matching*) a)
 304                            (setf *matching-path* (cddr *matching-path*))))
 305                     ))) ))
 306
 307
 308     (let ((matching ()))
 309       (maphash #'(lambda (u v) (if (< u v) (setq matching (cons `((mlist simp) ,u ,v) matching))))
 310                *matching-matching*)
 311       (cons '(mlist simp) matching))))
 312
 313 ;; returns the path from b to a
 314 (defun find-augmenting-path (a b)
 315   (cond
 316     ((= a b)
 317      (push a *matching-path*))
 318
 319     ((eql (gethash b *matching-state*) 'A)
 320      (push b *matching-path*)
 321      (push (gethash b *matching-matching*) *matching-path*)
 322      (find-augmenting-path a (gethash (gethash b *matching-matching*) *matching-prev*)))
 323
 324     (t
 325      (push b *matching-path*)
 326      (find-augmenting-path-1 (gethash b *matching-matching*) (first (gethash b *matching-bridge*)))
 327      (find-augmenting-path a (second (gethash b *matching-bridge*)))) ))
 328
 329 (defun find-augmenting-path-1 (a b)
 330   (cond
 331     ((= a b)
 332      (push a *matching-path*))
 333
 334     ((eql (gethash b *matching-state*) 'A)
 335      (find-augmenting-path-1 a (gethash (gethash b *matching-matching*) *matching-prev*))
 336      (push (gethash b *matching-matching*) *matching-path*)
 337      (push b *matching-path*))
 338
 339     (t
 340      (find-augmenting-path-1 a (second (gethash b *matching-bridge*)))
 341      (find-augmenting-path (gethash b *matching-matching*) (first (gethash b *matching-bridge*)))
 342      (push b *matching-path*)) ))