Index: .hgignore
===================================================================
--- .hgignore (revision 298)
+++ .hgignore (revision 400)
@@ -37,4 +37,5 @@
^objs.*/.*
^test/[a-z_]*$
+^tools/[a-z-_]*$
^demo/.*_demo$
^build/.*
Index: Makefile.am
===================================================================
--- Makefile.am (revision 321)
+++ Makefile.am (revision 375)
@@ -1,3 +1,5 @@
ACLOCAL_AMFLAGS = -I m4
+
+AM_CXXFLAGS = $(WARNINGCXXFLAGS)
AM_CPPFLAGS = -I$(top_srcdir) -I$(top_builddir)
Index: configure.ac
===================================================================
--- configure.ac (revision 310)
+++ configure.ac (revision 375)
@@ -19,6 +19,4 @@
AC_CONFIG_SRCDIR([lemon/list_graph.h])
AC_CONFIG_HEADERS([config.h lemon/config.h])
-
-lx_cmdline_cxxflags_set=${CXXFLAGS+set}
dnl Do compilation tests using the C++ compiler.
@@ -47,7 +45,8 @@
dnl Set custom compiler flags when using g++.
-if test x"$lx_cmdline_cxxflags_set" != x"set" -a "$GXX" = yes -a "$ICC" = no; then
- CXXFLAGS="$CXXFLAGS -Wall -W -Wall -W -Wunused -Wformat=2 -Wctor-dtor-privacy -Wnon-virtual-dtor -Wno-char-subscripts -Wwrite-strings -Wno-char-subscripts -Wreturn-type -Wcast-qual -Wcast-align -Wsign-promo -Woverloaded-virtual -Woverloaded-virtual -ansi -fno-strict-aliasing -Wold-style-cast -Wno-unknown-pragmas"
+if test "$GXX" = yes -a "$ICC" = no; then
+ WARNINGCXXFLAGS="-Wall -W -Wall -W -Wunused -Wformat=2 -Wctor-dtor-privacy -Wnon-virtual-dtor -Wno-char-subscripts -Wwrite-strings -Wno-char-subscripts -Wreturn-type -Wcast-qual -Wcast-align -Wsign-promo -Woverloaded-virtual -ansi -fno-strict-aliasing -Wold-style-cast -Wno-unknown-pragmas"
fi
+AC_SUBST([WARNINGCXXFLAGS])
dnl Checks for libraries.
@@ -114,5 +113,5 @@
echo
echo C++ compiler.................. : $CXX
-echo C++ compiles flags............ : $CXXFLAGS
+echo C++ compiles flags............ : $WARNINGCXXFLAGS $CXXFLAGS
echo
#echo GLPK support.................. : $lx_glpk_found
Index: doc/CMakeLists.txt
===================================================================
--- doc/CMakeLists.txt (revision 225)
+++ doc/CMakeLists.txt (revision 347)
@@ -14,4 +14,5 @@
COMMAND rm -rf gen-images
COMMAND mkdir gen-images
+ COMMAND ${GHOSTSCRIPT_EXECUTABLE} -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4 -sDEVICE=pngalpha -r18 -sOutputFile=gen-images/grid_graph.png ${CMAKE_CURRENT_SOURCE_DIR}/images/grid_graph.eps
COMMAND ${GHOSTSCRIPT_EXECUTABLE} -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4 -sDEVICE=pngalpha -r18 -sOutputFile=gen-images/nodeshape_0.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_0.eps
COMMAND ${GHOSTSCRIPT_EXECUTABLE} -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4 -sDEVICE=pngalpha -r18 -sOutputFile=gen-images/nodeshape_1.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_1.eps
Index: doc/Doxyfile.in
===================================================================
--- doc/Doxyfile.in (revision 316)
+++ doc/Doxyfile.in (revision 379)
@@ -67,5 +67,5 @@
ENABLED_SECTIONS =
MAX_INITIALIZER_LINES = 5
-SHOW_USED_FILES = YES
+SHOW_USED_FILES = NO
SHOW_DIRECTORIES = YES
SHOW_FILES = YES
Index: doc/Makefile.am
===================================================================
--- doc/Makefile.am (revision 317)
+++ doc/Makefile.am (revision 349)
@@ -15,4 +15,5 @@
DOC_EPS_IMAGES18 = \
+ grid_graph.eps \
nodeshape_0.eps \
nodeshape_1.eps \
Index: doc/groups.dox
===================================================================
--- doc/groups.dox (revision 318)
+++ doc/groups.dox (revision 403)
@@ -465,5 +465,5 @@
/**
-@defgroup lemon_io LEMON Input-Output
+@defgroup lemon_io LEMON Graph Format
@ingroup io_group
\brief Reading and writing LEMON Graph Format.
@@ -480,4 +480,20 @@
This group describes general \c EPS drawing methods and special
graph exporting tools.
+*/
+
+/**
+@defgroup dimacs_group DIMACS format
+@ingroup io_group
+\brief Read and write files in DIMACS format
+
+Tools to read a digraph from or write it to a file in DIMACS format data.
+*/
+
+/**
+@defgroup nauty_group NAUTY Format
+@ingroup io_group
+\brief Read \e Nauty format
+
+Tool to read graphs from \e Nauty format data.
*/
Index: doc/images/grid_graph.eps
===================================================================
--- doc/images/grid_graph.eps (revision 347)
+++ doc/images/grid_graph.eps (revision 347)
@@ -0,0 +1,286 @@
+%!PS-Adobe-2.0 EPSF-2.0
+%%Title: Grid undirected graph
+%%Copyright: (C) 2006 LEMON Project
+%%Creator: LEMON, graphToEps()
+%%CreationDate: Fri Sep 29 11:55:56 2006
+%%BoundingBox: 0 0 985 1144
+%%EndComments
+/lb { setlinewidth setrgbcolor newpath moveto
+ 4 2 roll 1 index 1 index curveto stroke } bind def
+/l { setlinewidth setrgbcolor newpath moveto lineto stroke } bind def
+/c { newpath dup 3 index add 2 index moveto 0 360 arc closepath } bind def
+/sq { newpath 2 index 1 index add 2 index 2 index add moveto
+ 2 index 1 index sub 2 index 2 index add lineto
+ 2 index 1 index sub 2 index 2 index sub lineto
+ 2 index 1 index add 2 index 2 index sub lineto
+ closepath pop pop pop} bind def
+/di { newpath 2 index 1 index add 2 index moveto
+ 2 index 2 index 2 index add lineto
+ 2 index 1 index sub 2 index lineto
+ 2 index 2 index 2 index sub lineto
+ closepath pop pop pop} bind def
+/nc { 0 0 0 setrgbcolor 5 index 5 index 5 index c fill
+ setrgbcolor 1.1 div c fill
+ } bind def
+/arrl 1 def
+/arrw 0.3 def
+/lrl { 2 index mul exch 2 index mul exch rlineto pop} bind def
+/arr { setrgbcolor /y1 exch def /x1 exch def /dy exch def /dx exch def
+ /w exch def /len exch def
+ newpath x1 dy w 2 div mul add y1 dx w 2 div mul sub moveto
+ len w sub arrl sub dx dy lrl
+ arrw dy dx neg lrl
+ dx arrl w add mul dy w 2 div arrw add mul sub
+ dy arrl w add mul dx w 2 div arrw add mul add rlineto
+ dx arrl w add mul neg dy w 2 div arrw add mul sub
+ dy arrl w add mul neg dx w 2 div arrw add mul add rlineto
+ arrw dy dx neg lrl
+ len w sub arrl sub neg dx dy lrl
+ closepath fill } bind def
+/cshow { 2 index 2 index moveto dup stringwidth pop
+ neg 2 div fosi .35 mul neg rmoveto show pop pop} def
+
+gsave
+2 2 scale
+50 40 translate
+5.5000 5.5000 scale
+% 1.14018 1.14018 translate
+%Edges:
+gsave
+70 80 70 90 0 0 0 0.5000 l
+70 70 70 80 0 0 0 0.5000 l
+70 60 70 70 0 0 0 0.5000 l
+70 50 70 60 0 0 0 0.5000 l
+70 40 70 50 0 0 0 0.5000 l
+70 30 70 40 0 0 0 0.5000 l
+70 20 70 30 0 0 0 0.5000 l
+70 10 70 20 0 0 0 0.5000 l
+70 0 70 10 0 0 0 0.5000 l
+60 80 60 90 0 0 0 0.5000 l
+60 70 60 80 0 0 0 0.5000 l
+60 60 60 70 0 0 0 0.5000 l
+60 50 60 60 0 0 0 0.5000 l
+60 40 60 50 0 0 0 0.5000 l
+60 30 60 40 0 0 0 0.5000 l
+60 20 60 30 0 0 0 0.5000 l
+60 10 60 20 0 0 0 0.5000 l
+60 0 60 10 0 0 0 0.5000 l
+50 80 50 90 0 0 0 0.5000 l
+50 70 50 80 0 0 0 0.5000 l
+50 60 50 70 0 0 0 0.5000 l
+50 50 50 60 0 0 0 0.5000 l
+50 40 50 50 0 0 0 0.5000 l
+50 30 50 40 0 0 0 0.5000 l
+50 20 50 30 0 0 0 0.5000 l
+50 10 50 20 0 0 0 0.5000 l
+50 0 50 10 0 0 0 0.5000 l
+40 80 40 90 0 0 0 0.5000 l
+40 70 40 80 0 0 0 0.5000 l
+40 60 40 70 0 0 0 0.5000 l
+40 50 40 60 0 0 0 0.5000 l
+40 40 40 50 0 0 0 0.5000 l
+40 30 40 40 0 0 0 0.5000 l
+40 20 40 30 0 0 0 0.5000 l
+40 10 40 20 0 0 0 0.5000 l
+40 0 40 10 0 0 0 0.5000 l
+30 80 30 90 0 0 0 0.5000 l
+30 70 30 80 0 0 0 0.5000 l
+30 60 30 70 0 0 0 0.5000 l
+30 50 30 60 0 0 0 0.5000 l
+30 40 30 50 0 0 0 0.5000 l
+30 30 30 40 0 0 0 0.5000 l
+30 20 30 30 0 0 0 0.5000 l
+30 10 30 20 0 0 0 0.5000 l
+30 0 30 10 0 0 0 0.5000 l
+20 80 20 90 0 0 0 0.5000 l
+20 70 20 80 0 0 0 0.5000 l
+20 60 20 70 0 0 0 0.5000 l
+20 50 20 60 0 0 0 0.5000 l
+20 40 20 50 0 0 0 0.5000 l
+20 30 20 40 0 0 0 0.5000 l
+20 20 20 30 0 0 0 0.5000 l
+20 10 20 20 0 0 0 0.5000 l
+20 0 20 10 0 0 0 0.5000 l
+10 80 10 90 0 0 0 0.5000 l
+10 70 10 80 0 0 0 0.5000 l
+10 60 10 70 0 0 0 0.5000 l
+10 50 10 60 0 0 0 0.5000 l
+10 40 10 50 0 0 0 0.5000 l
+10 30 10 40 0 0 0 0.5000 l
+10 20 10 30 0 0 0 0.5000 l
+10 10 10 20 0 0 0 0.5000 l
+10 0 10 10 0 0 0 0.5000 l
+0 80 0 90 0 0 0 0.5000 l
+0 70 0 80 0 0 0 0.5000 l
+0 60 0 70 0 0 0 0.5000 l
+0 50 0 60 0 0 0 0.5000 l
+0 40 0 50 0 0 0 0.5000 l
+0 30 0 40 0 0 0 0.5000 l
+0 20 0 30 0 0 0 0.5000 l
+0 10 0 20 0 0 0 0.5000 l
+0 0 0 10 0 0 0 0.5000 l
+60 90 70 90 0 0 0 0.5000 l
+60 80 70 80 0 0 0 0.5000 l
+60 70 70 70 0 0 0 0.5000 l
+60 60 70 60 0 0 0 0.5000 l
+60 50 70 50 0 0 0 0.5000 l
+60 40 70 40 0 0 0 0.5000 l
+60 30 70 30 0 0 0 0.5000 l
+60 20 70 20 0 0 0 0.5000 l
+60 10 70 10 0 0 0 0.5000 l
+60 0 70 0 0 0 0 0.5000 l
+50 90 60 90 0 0 0 0.5000 l
+50 80 60 80 0 0 0 0.5000 l
+50 70 60 70 0 0 0 0.5000 l
+50 60 60 60 0 0 0 0.5000 l
+50 50 60 50 0 0 0 0.5000 l
+50 40 60 40 0 0 0 0.5000 l
+50 30 60 30 0 0 0 0.5000 l
+50 20 60 20 0 0 0 0.5000 l
+50 10 60 10 0 0 0 0.5000 l
+50 0 60 0 0 0 0 0.5000 l
+40 90 50 90 0 0 0 0.5000 l
+40 80 50 80 0 0 0 0.5000 l
+40 70 50 70 0 0 0 0.5000 l
+40 60 50 60 0 0 0 0.5000 l
+40 50 50 50 0 0 0 0.5000 l
+40 40 50 40 0 0 0 0.5000 l
+40 30 50 30 0 0 0 0.5000 l
+40 20 50 20 0 0 0 0.5000 l
+40 10 50 10 0 0 0 0.5000 l
+40 0 50 0 0 0 0 0.5000 l
+30 90 40 90 0 0 0 0.5000 l
+30 80 40 80 0 0 0 0.5000 l
+30 70 40 70 0 0 0 0.5000 l
+30 60 40 60 0 0 0 0.5000 l
+30 50 40 50 0 0 0 0.5000 l
+30 40 40 40 0 0 0 0.5000 l
+30 30 40 30 0 0 0 0.5000 l
+30 20 40 20 0 0 0 0.5000 l
+30 10 40 10 0 0 0 0.5000 l
+30 0 40 0 0 0 0 0.5000 l
+20 90 30 90 0 0 0 0.5000 l
+20 80 30 80 0 0 0 0.5000 l
+20 70 30 70 0 0 0 0.5000 l
+20 60 30 60 0 0 0 0.5000 l
+20 50 30 50 0 0 0 0.5000 l
+20 40 30 40 0 0 0 0.5000 l
+20 30 30 30 0 0 0 0.5000 l
+20 20 30 20 0 0 0 0.5000 l
+20 10 30 10 0 0 0 0.5000 l
+20 0 30 0 0 0 0 0.5000 l
+10 90 20 90 0 0 0 0.5000 l
+10 80 20 80 0 0 0 0.5000 l
+10 70 20 70 0 0 0 0.5000 l
+10 60 20 60 0 0 0 0.5000 l
+10 50 20 50 0 0 0 0.5000 l
+10 40 20 40 0 0 0 0.5000 l
+10 30 20 30 0 0 0 0.5000 l
+10 20 20 20 0 0 0 0.5000 l
+10 10 20 10 0 0 0 0.5000 l
+10 0 20 0 0 0 0 0.5000 l
+0 90 10 90 0 0 0 0.5000 l
+0 80 10 80 0 0 0 0.5000 l
+0 70 10 70 0 0 0 0.5000 l
+0 60 10 60 0 0 0 0.5000 l
+0 50 10 50 0 0 0 0.5000 l
+0 40 10 40 0 0 0 0.5000 l
+0 30 10 30 0 0 0 0.5000 l
+0 20 10 20 0 0 0 0.5000 l
+0 10 10 10 0 0 0 0.5000 l
+0 0 10 0 0 0 0 0.5000 l
+grestore
+%Nodes:
+gsave
+70 90 1.4000 0 0 0 nc
+70 80 1.4000 1 1 1 nc
+70 70 1.4000 1 1 1 nc
+70 60 1.4000 1 1 1 nc
+70 50 1.4000 1 1 1 nc
+70 40 1.4000 1 1 1 nc
+70 30 1.4000 1 1 1 nc
+70 20 1.4000 1 1 1 nc
+70 10 1.4000 1 1 1 nc
+70 0 1.4000 0 0 0 nc
+60 90 1.4000 1 1 1 nc
+60 80 1.4000 1 1 1 nc
+60 70 1.4000 1 1 1 nc
+60 60 1.4000 1 1 1 nc
+60 50 1.4000 1 1 1 nc
+60 40 1.4000 1 1 1 nc
+60 30 1.4000 1 1 1 nc
+60 20 1.4000 1 1 1 nc
+60 10 1.4000 1 1 1 nc
+60 0 1.4000 1 1 1 nc
+50 90 1.4000 1 1 1 nc
+50 80 1.4000 1 1 1 nc
+50 70 1.4000 1 1 1 nc
+50 60 1.4000 1 1 1 nc
+50 50 1.4000 1 1 1 nc
+50 40 1.4000 1 1 1 nc
+50 30 1.4000 1 1 1 nc
+50 20 1.4000 1 1 1 nc
+50 10 1.4000 1 1 1 nc
+50 0 1.4000 1 1 1 nc
+40 90 1.4000 1 1 1 nc
+40 80 1.4000 1 1 1 nc
+40 70 1.4000 1 1 1 nc
+40 60 1.4000 1 1 1 nc
+40 50 1.4000 1 1 1 nc
+40 40 1.4000 1 1 1 nc
+40 30 1.4000 1 1 1 nc
+40 20 1.4000 1 1 1 nc
+40 10 1.4000 1 1 1 nc
+40 0 1.4000 1 1 1 nc
+30 90 1.4000 1 1 1 nc
+30 80 1.4000 1 1 1 nc
+30 70 1.4000 1 1 1 nc
+30 60 1.4000 1 1 1 nc
+30 50 1.4000 1 1 1 nc
+30 40 1.4000 1 1 1 nc
+30 30 1.4000 1 1 1 nc
+30 20 1.4000 1 1 1 nc
+30 10 1.4000 1 1 1 nc
+30 0 1.4000 1 1 1 nc
+20 90 1.4000 1 1 1 nc
+20 80 1.4000 1 1 1 nc
+20 70 1.4000 1 1 1 nc
+20 60 1.4000 1 1 1 nc
+20 50 1.4000 1 1 1 nc
+20 40 1.4000 1 1 1 nc
+20 30 1.4000 1 1 1 nc
+20 20 1.4000 1 1 1 nc
+20 10 1.4000 1 1 1 nc
+20 0 1.4000 1 1 1 nc
+10 90 1.4000 1 1 1 nc
+10 80 1.4000 1 1 1 nc
+10 70 1.4000 1 1 1 nc
+10 60 1.4000 1 1 1 nc
+10 50 1.4000 1 1 1 nc
+10 40 1.4000 1 1 1 nc
+10 30 1.4000 1 1 1 nc
+10 20 1.4000 1 1 1 nc
+10 10 1.4000 1 1 1 nc
+10 0 1.4000 1 1 1 nc
+0 90 1.4000 0 0 0 nc
+0 80 1.4000 1 1 1 nc
+0 70 1.4000 1 1 1 nc
+0 60 1.4000 1 1 1 nc
+0 50 1.4000 1 1 1 nc
+0 40 1.4000 1 1 1 nc
+0 30 1.4000 1 1 1 nc
+0 20 1.4000 1 1 1 nc
+0 10 1.4000 1 1 1 nc
+0 0 1.4000 0 0 0 nc
+grestore
+gsave
+/fosi 3.5 def
+(Helvetica) findfont fosi scalefont setfont
+0 0 0 setrgbcolor
+0 95 ((0,height-1)) cshow
+67 95 ((width-1,height-1)) cshow
+0 -5 ((0,0)) cshow
+70 -5 ((width-1,0)) cshow
+grestore
+grestore
+showpage
Index: doc/migration.dox
===================================================================
--- doc/migration.dox (revision 314)
+++ doc/migration.dox (revision 355)
@@ -26,5 +26,5 @@
Many of these changes adjusted automatically by the
-script/lemon-0.x-to-1.x.sh tool. Those requiring manual
+lemon-0.x-to-1.x.sh tool. Those requiring manual
update are typeset boldface.
@@ -54,7 +54,9 @@
\warning
-The script/lemon-0.x-to-1.x.sh tool replaces all instances of
-the words \c graph, \c digraph, \c edge and \c arc, so it replaces them
-in strings, comments etc. as well as in all identifiers.
+The lemon-0.x-to-1.x.sh script replaces the words \c graph,
+\c ugraph, \c edge and \c uedge in your own identifiers and in
+strings, comments etc. as well as in all LEMON specific identifiers.
+So use the script carefully and make a backup copy of your source files
+before applying the script to them.
\section migration-lgf LGF tools
Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am (revision 259)
+++ lemon/Makefile.am (revision 410)
@@ -13,5 +13,5 @@
lemon/random.cc
-#lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS)
+#lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS) $(AM_CXXFLAGS)
#lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS)
@@ -28,6 +28,11 @@
lemon/dijkstra.h \
lemon/dim2.h \
+ lemon/dimacs.h \
+ lemon/elevator.h \
lemon/error.h \
+ lemon/full_graph.h \
lemon/graph_to_eps.h \
+ lemon/grid_graph.h \
+ lemon/hypercube_graph.h \
lemon/kruskal.h \
lemon/lgf_reader.h \
@@ -36,7 +41,11 @@
lemon/maps.h \
lemon/math.h \
+ lemon/max_matching.h \
+ lemon/nauty_reader.h \
lemon/path.h \
+ lemon/preflow.h \
lemon/random.h \
lemon/smart_graph.h \
+ lemon/suurballe.h \
lemon/time_measure.h \
lemon/tolerance.h \
Index: lemon/bfs.h
===================================================================
--- lemon/bfs.h (revision 301)
+++ lemon/bfs.h (revision 341)
@@ -52,5 +52,5 @@
///Instantiates a PredMap.
- ///This function instantiates a PredMap.
+ ///This function instantiates a PredMap.
///\param g is the digraph, to which we would like to define the
///PredMap.
@@ -81,6 +81,5 @@
///The type of the map that indicates which nodes are reached.
- ///The type of the map that indicates which nodes are reached.
- ///It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+ ///The type of the map that indicates which nodes are reached.///It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
typedef typename Digraph::template NodeMap ReachedMap;
///Instantiates a ReachedMap.
Index: lemon/bits/alteration_notifier.h
===================================================================
--- lemon/bits/alteration_notifier.h (revision 314)
+++ lemon/bits/alteration_notifier.h (revision 373)
@@ -36,59 +36,60 @@
// a container.
//
- // The simple graph's can be refered as two containers, one node container
- // and one edge container. But they are not standard containers they
- // does not store values directly they are just key continars for more
- // value containers which are the node and edge maps.
- //
- // The graph's node and edge sets can be changed as we add or erase
+ // The simple graphs can be refered as two containers: a node container
+ // and an edge container. But they do not store values directly, they
+ // are just key continars for more value containers, which are the
+ // node and edge maps.
+ //
+ // The node and edge sets of the graphs can be changed as we add or erase
// nodes and edges in the graph. LEMON would like to handle easily
// that the node and edge maps should contain values for all nodes or
// edges. If we want to check on every indicing if the map contains
// the current indicing key that cause a drawback in the performance
- // in the library. We use another solution we notify all maps about
+ // in the library. We use another solution: we notify all maps about
// an alteration in the graph, which cause only drawback on the
// alteration of the graph.
//
- // This class provides an interface to the container. The \e first() and \e
- // next() member functions make possible to iterate on the keys of the
- // container. The \e id() function returns an integer id for each key.
- // The \e maxId() function gives back an upper bound of the ids.
+ // This class provides an interface to a node or edge container.
+ // The first() and next() member functions make possible
+ // to iterate on the keys of the container.
+ // The id() function returns an integer id for each key.
+ // The maxId() function gives back an upper bound of the ids.
//
// For the proper functonality of this class, we should notify it
- // about each alteration in the container. The alterations have four type
- // as \e add(), \e erase(), \e build() and \e clear(). The \e add() and
- // \e erase() signals that only one or few items added or erased to or
- // from the graph. If all items are erased from the graph or from an empty
- // graph a new graph is builded then it can be signaled with the
+ // about each alteration in the container. The alterations have four type:
+ // add(), erase(), build() and clear(). The add() and
+ // erase() signal that only one or few items added or erased to or
+ // from the graph. If all items are erased from the graph or if a new graph
+ // is built from an empty graph, then it can be signaled with the
// clear() and build() members. Important rule that if we erase items
- // from graph we should first signal the alteration and after that erase
+ // from graphs we should first signal the alteration and after that erase
// them from the container, on the other way on item addition we should
// first extend the container and just after that signal the alteration.
//
// The alteration can be observed with a class inherited from the
- // \e ObserverBase nested class. The signals can be handled with
+ // ObserverBase nested class. The signals can be handled with
// overriding the virtual functions defined in the base class. The
// observer base can be attached to the notifier with the
- // \e attach() member and can be detached with detach() function. The
+ // attach() member and can be detached with detach() function. The
// alteration handlers should not call any function which signals
// an other alteration in the same notifier and should not
// detach any observer from the notifier.
//
- // Alteration observers try to be exception safe. If an \e add() or
- // a \e clear() function throws an exception then the remaining
+ // Alteration observers try to be exception safe. If an add() or
+ // a clear() function throws an exception then the remaining
// observeres will not be notified and the fulfilled additions will
- // be rolled back by calling the \e erase() or \e clear()
- // functions. Thence the \e erase() and \e clear() should not throw
- // exception. Actullay, it can be throw only \ref ImmediateDetach
- // exception which detach the observer from the notifier.
- //
- // There are some place when the alteration observing is not completly
+ // be rolled back by calling the erase() or clear() functions.
+ // Hence erase() and clear() should not throw exception.
+ // Actullay, they can throw only \ref ImmediateDetach exception,
+ // which detach the observer from the notifier.
+ //
+ // There are some cases, when the alteration observing is not completly
// reliable. If we want to carry out the node degree in the graph
- // as in the \ref InDegMap and we use the reverseEdge that cause
+ // as in the \ref InDegMap and we use the reverseArc(), then it cause
// unreliable functionality. Because the alteration observing signals
- // only erasing and adding but not the reversing it will stores bad
- // degrees. The sub graph adaptors cannot signal the alterations because
- // just a setting in the filter map can modify the graph and this cannot
- // be watched in any way.
+ // only erasing and adding but not the reversing, it will stores bad
+ // degrees. Apart form that the subgraph adaptors cannot even signal
+ // the alterations because just a setting in the filter map can modify
+ // the graph and this cannot be watched in any way.
//
// \param _Container The container which is observed.
@@ -104,11 +105,11 @@
typedef _Item Item;
- // \brief Exception which can be called from \e clear() and
- // \e erase().
- //
- // From the \e clear() and \e erase() function only this
+ // \brief Exception which can be called from clear() and
+ // erase().
+ //
+ // From the clear() and erase() function only this
// exception is allowed to throw. The exception immediatly
// detaches the current observer from the notifier. Because the
- // \e clear() and \e erase() should not throw other exceptions
+ // clear() and erase() should not throw other exceptions
// it can be used to invalidate the observer.
struct ImmediateDetach {};
@@ -122,6 +123,5 @@
// The observer interface contains some pure virtual functions
// to override. The add() and erase() functions are
- // to notify the oberver when one item is added or
- // erased.
+ // to notify the oberver when one item is added or erased.
//
// The build() and clear() members are to notify the observer
Index: lemon/bits/array_map.h
===================================================================
--- lemon/bits/array_map.h (revision 314)
+++ lemon/bits/array_map.h (revision 373)
@@ -37,10 +37,9 @@
// \brief Graph map based on the array storage.
//
- // The ArrayMap template class is graph map structure what
- // automatically updates the map when a key is added to or erased from
- // the map. This map uses the allocators to implement
- // the container functionality.
+ // The ArrayMap template class is graph map structure that automatically
+ // updates the map when a key is added to or erased from the graph.
+ // This map uses the allocators to implement the container functionality.
//
- // The template parameters are the Graph the current Item type and
+ // The template parameters are the Graph, the current Item type and
// the Value type of the map.
template
@@ -48,12 +47,12 @@
: public ItemSetTraits<_Graph, _Item>::ItemNotifier::ObserverBase {
public:
- // The graph type of the maps.
+ // The graph type.
typedef _Graph Graph;
- // The item type of the map.
+ // The item type.
typedef _Item Item;
// The reference map tag.
typedef True ReferenceMapTag;
- // The key type of the maps.
+ // The key type of the map.
typedef _Item Key;
// The value type of the map.
@@ -201,5 +200,5 @@
// \brief Adds a new key to the map.
//
- // It adds a new key to the map. It called by the observer notifier
+ // It adds a new key to the map. It is called by the observer notifier
// and it overrides the add() member function of the observer base.
virtual void add(const Key& key) {
@@ -229,5 +228,5 @@
// \brief Adds more new keys to the map.
//
- // It adds more new keys to the map. It called by the observer notifier
+ // It adds more new keys to the map. It is called by the observer notifier
// and it overrides the add() member function of the observer base.
virtual void add(const std::vector& keys) {
@@ -273,5 +272,5 @@
// \brief Erase a key from the map.
//
- // Erase a key from the map. It called by the observer notifier
+ // Erase a key from the map. It is called by the observer notifier
// and it overrides the erase() member function of the observer base.
virtual void erase(const Key& key) {
@@ -282,5 +281,5 @@
// \brief Erase more keys from the map.
//
- // Erase more keys from the map. It called by the observer notifier
+ // Erase more keys from the map. It is called by the observer notifier
// and it overrides the erase() member function of the observer base.
virtual void erase(const std::vector& keys) {
@@ -291,7 +290,7 @@
}
- // \brief Buildes the map.
- //
- // It buildes the map. It called by the observer notifier
+ // \brief Builds the map.
+ //
+ // It builds the map. It is called by the observer notifier
// and it overrides the build() member function of the observer base.
virtual void build() {
@@ -307,5 +306,5 @@
// \brief Clear the map.
//
- // It erase all items from the map. It called by the observer notifier
+ // It erase all items from the map. It is called by the observer notifier
// and it overrides the clear() member function of the observer base.
virtual void clear() {
Index: lemon/bits/base_extender.h
===================================================================
--- lemon/bits/base_extender.h (revision 314)
+++ lemon/bits/base_extender.h (revision 373)
@@ -31,5 +31,5 @@
//\ingroup digraphbits
//\file
-//\brief Extenders for the digraph types
+//\brief Extenders for the graph types
namespace lemon {
Index: lemon/bits/graph_extender.h
===================================================================
--- lemon/bits/graph_extender.h (revision 314)
+++ lemon/bits/graph_extender.h (revision 373)
@@ -30,10 +30,10 @@
//\ingroup graphbits
//\file
-//\brief Extenders for the digraph types
+//\brief Extenders for the graph types
namespace lemon {
// \ingroup graphbits
//
- // \brief Extender for the Digraphs
+ // \brief Extender for the digraph implementations
template
class DigraphExtender : public Base {
Index: lemon/bits/traits.h
===================================================================
--- lemon/bits/traits.h (revision 314)
+++ lemon/bits/traits.h (revision 372)
@@ -219,4 +219,17 @@
template
+ struct ArcNumTagIndicator {
+ static const bool value = false;
+ };
+
+ template
+ struct ArcNumTagIndicator<
+ Graph,
+ typename enable_if::type
+ > {
+ static const bool value = true;
+ };
+
+ template
struct EdgeNumTagIndicator {
static const bool value = false;
@@ -232,4 +245,17 @@
template
+ struct FindArcTagIndicator {
+ static const bool value = false;
+ };
+
+ template
+ struct FindArcTagIndicator<
+ Graph,
+ typename enable_if::type
+ > {
+ static const bool value = true;
+ };
+
+ template
struct FindEdgeTagIndicator {
static const bool value = false;
Index: lemon/bits/vector_map.h
===================================================================
--- lemon/bits/vector_map.h (revision 314)
+++ lemon/bits/vector_map.h (revision 373)
@@ -39,7 +39,7 @@
// \brief Graph map based on the std::vector storage.
//
- // The VectorMap template class is graph map structure what
- // automatically updates the map when a key is added to or erased from
- // the map. This map type uses the std::vector to store the values.
+ // The VectorMap template class is graph map structure that automatically
+ // updates the map when a key is added to or erased from the graph.
+ // This map type uses std::vector to store the values.
//
// \tparam _Graph The graph this map is attached to.
@@ -170,5 +170,5 @@
// \brief Adds a new key to the map.
//
- // It adds a new key to the map. It called by the observer notifier
+ // It adds a new key to the map. It is called by the observer notifier
// and it overrides the add() member function of the observer base.
virtual void add(const Key& key) {
@@ -181,5 +181,5 @@
// \brief Adds more new keys to the map.
//
- // It adds more new keys to the map. It called by the observer notifier
+ // It adds more new keys to the map. It is called by the observer notifier
// and it overrides the add() member function of the observer base.
virtual void add(const std::vector& keys) {
@@ -196,5 +196,5 @@
// \brief Erase a key from the map.
//
- // Erase a key from the map. It called by the observer notifier
+ // Erase a key from the map. It is called by the observer notifier
// and it overrides the erase() member function of the observer base.
virtual void erase(const Key& key) {
@@ -204,5 +204,5 @@
// \brief Erase more keys from the map.
//
- // Erase more keys from the map. It called by the observer notifier
+ // It erases more keys from the map. It is called by the observer notifier
// and it overrides the erase() member function of the observer base.
virtual void erase(const std::vector& keys) {
@@ -212,7 +212,7 @@
}
- // \brief Buildes the map.
- //
- // It buildes the map. It called by the observer notifier
+ // \brief Build the map.
+ //
+ // It builds the map. It is called by the observer notifier
// and it overrides the build() member function of the observer base.
virtual void build() {
@@ -224,5 +224,5 @@
// \brief Clear the map.
//
- // It erase all items from the map. It called by the observer notifier
+ // It erases all items from the map. It is called by the observer notifier
// and it overrides the clear() member function of the observer base.
virtual void clear() {
Index: lemon/dimacs.h
===================================================================
--- lemon/dimacs.h (revision 403)
+++ lemon/dimacs.h (revision 403)
@@ -0,0 +1,357 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_DIMACS_H
+#define LEMON_DIMACS_H
+
+#include
+#include
+#include
+#include
+#include
+
+/// \ingroup dimacs_group
+/// \file
+/// \brief DIMACS file format reader.
+
+namespace lemon {
+
+ /// \addtogroup dimacs_group
+ /// @{
+
+ /// DIMACS file type descriptor.
+ struct DimacsDescriptor
+ {
+ ///File type enum
+ enum Type
+ {
+ NONE, MIN, MAX, SP, MAT
+ };
+ ///The file type
+ Type type;
+ ///The number of nodes in the graph
+ int nodeNum;
+ ///The number of edges in the graph
+ int edgeNum;
+ int lineShift;
+ /// Constructor. Sets the type to NONE.
+ DimacsDescriptor() : type(NONE) {}
+ };
+
+ ///Discover the type of a DIMACS file
+
+ ///It starts seeking the begining of the file for the problem type
+ ///and size info. The found data is returned in a special struct
+ ///that can be evaluated and passed to the appropriate reader
+ ///function.
+ DimacsDescriptor dimacsType(std::istream& is)
+ {
+ DimacsDescriptor r;
+ std::string problem,str;
+ char c;
+ r.lineShift=0;
+ while (is >> c)
+ switch(c)
+ {
+ case 'p':
+ if(is >> problem >> r.nodeNum >> r.edgeNum)
+ {
+ getline(is, str);
+ r.lineShift++;
+ if(problem=="min") r.type=DimacsDescriptor::MIN;
+ else if(problem=="max") r.type=DimacsDescriptor::MAX;
+ else if(problem=="sp") r.type=DimacsDescriptor::SP;
+ else if(problem=="mat") r.type=DimacsDescriptor::MAT;
+ else throw FormatError("Unknown problem type");
+ return r;
+ }
+ else
+ {
+ throw FormatError("Missing or wrong problem type declaration.");
+ }
+ break;
+ case 'c':
+ getline(is, str);
+ r.lineShift++;
+ break;
+ default:
+ throw FormatError("Unknown DIMACS declaration.");
+ }
+ throw FormatError("Missing problem type declaration.");
+ }
+
+
+
+ /// DIMACS minimum cost flow reader function.
+ ///
+ /// This function reads a minimum cost flow instance from DIMACS format,
+ /// i.e. from a DIMACS file having a line starting with
+ /// \code
+ /// p min
+ /// \endcode
+ /// At the beginning, \c g is cleared by \c g.clear(). The supply
+ /// amount of the nodes are written to \c supply (signed). The
+ /// lower bounds, capacities and costs of the arcs are written to
+ /// \c lower, \c capacity and \c cost.
+ ///
+ /// If the file type was previously evaluated by dimacsType(), then
+ /// the descriptor struct should be given by the \c dest parameter.
+ template
+ void readDimacsMin(std::istream& is,
+ Digraph &g,
+ LowerMap& lower,
+ CapacityMap& capacity,
+ CostMap& cost,
+ SupplyMap& supply,
+ DimacsDescriptor desc=DimacsDescriptor())
+ {
+ g.clear();
+ std::vector nodes;
+ typename Digraph::Arc e;
+ std::string problem, str;
+ char c;
+ int i, j;
+ if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+ if(desc.type!=DimacsDescriptor::MIN)
+ throw FormatError("Problem type mismatch");
+
+ nodes.resize(desc.nodeNum + 1);
+ for (int k = 1; k <= desc.nodeNum; ++k) {
+ nodes[k] = g.addNode();
+ supply.set(nodes[k], 0);
+ }
+
+ typename SupplyMap::Value sup;
+ typename CapacityMap::Value low;
+ typename CapacityMap::Value cap;
+ typename CostMap::Value co;
+ while (is >> c) {
+ switch (c) {
+ case 'c': // comment line
+ getline(is, str);
+ break;
+ case 'n': // node definition line
+ is >> i >> sup;
+ getline(is, str);
+ supply.set(nodes[i], sup);
+ break;
+ case 'a': // arc (arc) definition line
+ is >> i >> j >> low >> cap >> co;
+ getline(is, str);
+ e = g.addArc(nodes[i], nodes[j]);
+ lower.set(e, low);
+ if (cap >= 0)
+ capacity.set(e, cap);
+ else
+ capacity.set(e, -1);
+ cost.set(e, co);
+ break;
+ }
+ }
+ }
+
+ template
+ void _readDimacs(std::istream& is,
+ Digraph &g,
+ CapacityMap& capacity,
+ typename Digraph::Node &s,
+ typename Digraph::Node &t,
+ DimacsDescriptor desc=DimacsDescriptor()) {
+ g.clear();
+ s=t=INVALID;
+ std::vector nodes;
+ typename Digraph::Arc e;
+ char c, d;
+ int i, j;
+ typename CapacityMap::Value _cap;
+ std::string str;
+ nodes.resize(desc.nodeNum + 1);
+ for (int k = 1; k <= desc.nodeNum; ++k) {
+ nodes[k] = g.addNode();
+ }
+
+ while (is >> c) {
+ switch (c) {
+ case 'c': // comment line
+ getline(is, str);
+ break;
+ case 'n': // node definition line
+ if (desc.type==DimacsDescriptor::SP) { // shortest path problem
+ is >> i;
+ getline(is, str);
+ s = nodes[i];
+ }
+ if (desc.type==DimacsDescriptor::MAX) { // max flow problem
+ is >> i >> d;
+ getline(is, str);
+ if (d == 's') s = nodes[i];
+ if (d == 't') t = nodes[i];
+ }
+ break;
+ case 'a': // arc (arc) definition line
+ if (desc.type==DimacsDescriptor::SP ||
+ desc.type==DimacsDescriptor::MAX) {
+ is >> i >> j >> _cap;
+ getline(is, str);
+ e = g.addArc(nodes[i], nodes[j]);
+ capacity.set(e, _cap);
+ } else {
+ is >> i >> j;
+ getline(is, str);
+ g.addArc(nodes[i], nodes[j]);
+ }
+ break;
+ }
+ }
+ }
+
+ /// DIMACS maximum flow reader function.
+ ///
+ /// This function reads a maximum flow instance from DIMACS format,
+ /// i.e. from a DIMACS file having a line starting with
+ /// \code
+ /// p max
+ /// \endcode
+ /// At the beginning, \c g is cleared by \c g.clear(). The arc
+ /// capacities are written to \c capacity and \c s and \c t are
+ /// set to the source and the target nodes.
+ ///
+ /// If the file type was previously evaluated by dimacsType(), then
+ /// the descriptor struct should be given by the \c dest parameter.
+ template
+ void readDimacsMax(std::istream& is,
+ Digraph &g,
+ CapacityMap& capacity,
+ typename Digraph::Node &s,
+ typename Digraph::Node &t,
+ DimacsDescriptor desc=DimacsDescriptor()) {
+ if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+ if(desc.type!=DimacsDescriptor::MAX)
+ throw FormatError("Problem type mismatch");
+ _readDimacs(is,g,capacity,s,t,desc);
+ }
+
+ /// DIMACS shortest path reader function.
+ ///
+ /// This function reads a shortest path instance from DIMACS format,
+ /// i.e. from a DIMACS file having a line starting with
+ /// \code
+ /// p sp
+ /// \endcode
+ /// At the beginning, \c g is cleared by \c g.clear(). The arc
+ /// lengths are written to \c length and \c s is set to the
+ /// source node.
+ ///
+ /// If the file type was previously evaluated by dimacsType(), then
+ /// the descriptor struct should be given by the \c dest parameter.
+ template
+ void readDimacsSp(std::istream& is,
+ Digraph &g,
+ LengthMap& length,
+ typename Digraph::Node &s,
+ DimacsDescriptor desc=DimacsDescriptor()) {
+ typename Digraph::Node t;
+ if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+ if(desc.type!=DimacsDescriptor::SP)
+ throw FormatError("Problem type mismatch");
+ _readDimacs(is, g, length, s, t,desc);
+ }
+
+ /// DIMACS capacitated digraph reader function.
+ ///
+ /// This function reads an arc capacitated digraph instance from
+ /// DIMACS 'mat' or 'sp' format.
+ /// At the beginning, \c g is cleared by \c g.clear()
+ /// and the arc capacities/lengths are written to \c capacity.
+ ///
+ /// If the file type was previously evaluated by dimacsType(), then
+ /// the descriptor struct should be given by the \c dest parameter.
+ template
+ void readDimacsCap(std::istream& is,
+ Digraph &g,
+ CapacityMap& capacity,
+ DimacsDescriptor desc=DimacsDescriptor()) {
+ typename Digraph::Node u,v;
+ if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+ if(desc.type!=DimacsDescriptor::MAX || desc.type!=DimacsDescriptor::SP)
+ throw FormatError("Problem type mismatch");
+ _readDimacs(is, g, capacity, u, v, desc);
+ }
+
+ /// DIMACS plain digraph reader function.
+ ///
+ /// This function reads a digraph without any designated nodes and
+ /// maps from DIMACS format, i.e. from DIMACS files having a line
+ /// starting with
+ /// \code
+ /// p mat
+ /// \endcode
+ /// At the beginning, \c g is cleared by \c g.clear().
+ ///
+ /// If the file type was previously evaluated by dimacsType(), then
+ /// the descriptor struct should be given by the \c dest parameter.
+ template
+ void readDimacsMat(std::istream& is, Digraph &g,
+ DimacsDescriptor desc=DimacsDescriptor()) {
+ typename Digraph::Node u,v;
+ NullMap n;
+ if(desc.type==DimacsDescriptor::NONE) desc=dimacsType(is);
+ if(desc.type!=DimacsDescriptor::MAT)
+ throw FormatError("Problem type mismatch");
+ _readDimacs(is, g, n, u, v, desc);
+ }
+
+ /// DIMACS plain digraph writer function.
+ ///
+ /// This function writes a digraph without any designated nodes and
+ /// maps into DIMACS format, i.e. into DIMACS file having a line
+ /// starting with
+ /// \code
+ /// p mat
+ /// \endcode
+ /// If \c comment is not empty, then it will be printed in the first line
+ /// prefixed by 'c'.
+ template
+ void writeDimacsMat(std::ostream& os, const Digraph &g,
+ std::string comment="") {
+ typedef typename Digraph::NodeIt NodeIt;
+ typedef typename Digraph::ArcIt ArcIt;
+
+ if(!comment.empty())
+ os << "c " << comment << std::endl;
+ os << "p mat " << g.nodeNum() << " " << g.arcNum() << std::endl;
+
+ typename Digraph::template NodeMap nodes(g);
+ int i = 1;
+ for(NodeIt v(g); v != INVALID; ++v) {
+ nodes.set(v, i);
+ ++i;
+ }
+ for(ArcIt e(g); e != INVALID; ++e) {
+ os << "a " << nodes[g.source(e)] << " " << nodes[g.target(e)]
+ << std::endl;
+ }
+ }
+
+ /// @}
+
+} //namespace lemon
+
+#endif //LEMON_DIMACS_H
Index: lemon/elevator.h
===================================================================
--- lemon/elevator.h (revision 398)
+++ lemon/elevator.h (revision 398)
@@ -0,0 +1,981 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_ELEVATOR_H
+#define LEMON_ELEVATOR_H
+
+///\ingroup auxdat
+///\file
+///\brief Elevator class
+///
+///Elevator class implements an efficient data structure
+///for labeling items in push-relabel type algorithms.
+///
+
+#include
+
+namespace lemon {
+
+ ///Class for handling "labels" in push-relabel type algorithms.
+
+ ///A class for handling "labels" in push-relabel type algorithms.
+ ///
+ ///\ingroup auxdat
+ ///Using this class you can assign "labels" (nonnegative integer numbers)
+ ///to the edges or nodes of a graph, manipulate and query them through
+ ///operations typically arising in "push-relabel" type algorithms.
+ ///
+ ///Each item is either \em active or not, and you can also choose a
+ ///highest level active item.
+ ///
+ ///\sa LinkedElevator
+ ///
+ ///\param Graph Type of the underlying graph.
+ ///\param Item Type of the items the data is assigned to (Graph::Node,
+ ///Graph::Arc, Graph::Edge).
+ template
+ class Elevator
+ {
+ public:
+
+ typedef Item Key;
+ typedef int Value;
+
+ private:
+
+ typedef Item *Vit;
+ typedef typename ItemSetTraits::template Map::Type VitMap;
+ typedef typename ItemSetTraits::template Map::Type IntMap;
+
+ const Graph &_g;
+ int _max_level;
+ int _item_num;
+ VitMap _where;
+ IntMap _level;
+ std::vector- _items;
+ std::vector _first;
+ std::vector _last_active;
+
+ int _highest_active;
+
+ void copy(Item i, Vit p)
+ {
+ _where.set(*p=i,p);
+ }
+ void copy(Vit s, Vit p)
+ {
+ if(s!=p)
+ {
+ Item i=*s;
+ *p=i;
+ _where.set(i,p);
+ }
+ }
+ void swap(Vit i, Vit j)
+ {
+ Item ti=*i;
+ Vit ct = _where[ti];
+ _where.set(ti,_where[*i=*j]);
+ _where.set(*j,ct);
+ *j=ti;
+ }
+
+ public:
+
+ ///Constructor with given maximum level.
+
+ ///Constructor with given maximum level.
+ ///
+ ///\param graph The underlying graph.
+ ///\param max_level The maximum allowed level.
+ ///Set the range of the possible labels to [0..max_level].
+ Elevator(const Graph &graph,int max_level) :
+ _g(graph),
+ _max_level(max_level),
+ _item_num(_max_level),
+ _where(graph),
+ _level(graph,0),
+ _items(_max_level),
+ _first(_max_level+2),
+ _last_active(_max_level+2),
+ _highest_active(-1) {}
+ ///Constructor.
+
+ ///Constructor.
+ ///
+ ///\param graph The underlying graph.
+ ///Set the range of the possible labels to [0..max_level],
+ ///where \c max_level is equal to the number of labeled items in the graph.
+ Elevator(const Graph &graph) :
+ _g(graph),
+ _max_level(countItems(graph)),
+ _item_num(_max_level),
+ _where(graph),
+ _level(graph,0),
+ _items(_max_level),
+ _first(_max_level+2),
+ _last_active(_max_level+2),
+ _highest_active(-1)
+ {
+ }
+
+ ///Activate item \c i.
+
+ ///Activate item \c i.
+ ///\pre Item \c i shouldn't be active before.
+ void activate(Item i)
+ {
+ const int l=_level[i];
+ swap(_where[i],++_last_active[l]);
+ if(l>_highest_active) _highest_active=l;
+ }
+
+ ///Deactivate item \c i.
+
+ ///Deactivate item \c i.
+ ///\pre Item \c i must be active before.
+ void deactivate(Item i)
+ {
+ swap(_where[i],_last_active[_level[i]]--);
+ while(_highest_active>=0 &&
+ _last_active[_highest_active]<_first[_highest_active])
+ _highest_active--;
+ }
+
+ ///Query whether item \c i is active
+ bool active(Item i) const { return _where[i]<=_last_active[_level[i]]; }
+
+ ///Return the level of item \c i.
+ int operator[](Item i) const { return _level[i]; }
+
+ ///Return the number of items on level \c l.
+ int onLevel(int l) const
+ {
+ return _first[l+1]-_first[l];
+ }
+ ///Return true if level \c l is empty.
+ bool emptyLevel(int l) const
+ {
+ return _first[l+1]-_first[l]==0;
+ }
+ ///Return the number of items above level \c l.
+ int aboveLevel(int l) const
+ {
+ return _first[_max_level+1]-_first[l+1];
+ }
+ ///Return the number of active items on level \c l.
+ int activesOnLevel(int l) const
+ {
+ return _last_active[l]-_first[l]+1;
+ }
+ ///Return true if there is no active item on level \c l.
+ bool activeFree(int l) const
+ {
+ return _last_active[l]<_first[l];
+ }
+ ///Return the maximum allowed level.
+ int maxLevel() const
+ {
+ return _max_level;
+ }
+
+ ///\name Highest Active Item
+ ///Functions for working with the highest level
+ ///active item.
+
+ ///@{
+
+ ///Return a highest level active item.
+
+ ///Return a highest level active item or INVALID if there is no active
+ ///item.
+ Item highestActive() const
+ {
+ return _highest_active>=0?*_last_active[_highest_active]:INVALID;
+ }
+
+ ///Return the highest active level.
+
+ ///Return the level of the highest active item or -1 if there is no active
+ ///item.
+ int highestActiveLevel() const
+ {
+ return _highest_active;
+ }
+
+ ///Lift the highest active item by one.
+
+ ///Lift the item returned by highestActive() by one.
+ ///
+ void liftHighestActive()
+ {
+ Item it = *_last_active[_highest_active];
+ _level.set(it,_level[it]+1);
+ swap(_last_active[_highest_active]--,_last_active[_highest_active+1]);
+ --_first[++_highest_active];
+ }
+
+ ///Lift the highest active item to the given level.
+
+ ///Lift the item returned by highestActive() to level \c new_level.
+ ///
+ ///\warning \c new_level must be strictly higher
+ ///than the current level.
+ ///
+ void liftHighestActive(int new_level)
+ {
+ const Item li = *_last_active[_highest_active];
+
+ copy(--_first[_highest_active+1],_last_active[_highest_active]--);
+ for(int l=_highest_active+1;l=0 &&
+ _last_active[_highest_active]<_first[_highest_active])
+ _highest_active--;
+ }
+
+ ///@}
+
+ ///\name Active Item on Certain Level
+ ///Functions for working with the active items.
+
+ ///@{
+
+ ///Return an active item on level \c l.
+
+ ///Return an active item on level \c l or \ref INVALID if there is no such
+ ///an item. (\c l must be from the range [0...\c max_level].
+ Item activeOn(int l) const
+ {
+ return _last_active[l]>=_first[l]?*_last_active[l]:INVALID;
+ }
+
+ ///Lift the active item returned by \c activeOn(level) by one.
+
+ ///Lift the active item returned by \ref activeOn() "activeOn(level)"
+ ///by one.
+ Item liftActiveOn(int level)
+ {
+ Item it =*_last_active[level];
+ _level.set(it,_level[it]+1);
+ swap(_last_active[level]--, --_first[level+1]);
+ if (level+1>_highest_active) ++_highest_active;
+ }
+
+ ///Lift the active item returned by \c activeOn(level) to the given level.
+
+ ///Lift the active item returned by \ref activeOn() "activeOn(level)"
+ ///to the given level.
+ void liftActiveOn(int level, int new_level)
+ {
+ const Item ai = *_last_active[level];
+
+ copy(--_first[level+1], _last_active[level]--);
+ for(int l=level+1;l_highest_active) _highest_active=new_level;
+ }
+
+ ///Lift the active item returned by \c activeOn(level) to the top level.
+
+ ///Lift the active item returned by \ref activeOn() "activeOn(level)"
+ ///to the top level and deactivate it.
+ void liftActiveToTop(int level)
+ {
+ const Item ai = *_last_active[level];
+
+ copy(--_first[level+1],_last_active[level]--);
+ for(int l=level+1;l<_max_level;l++)
+ {
+ copy(_last_active[l],_first[l]);
+ copy(--_first[l+1], _last_active[l]--);
+ }
+ copy(ai,_first[_max_level]);
+ --_last_active[_max_level];
+ _level.set(ai,_max_level);
+
+ if (_highest_active==level) {
+ while(_highest_active>=0 &&
+ _last_active[_highest_active]<_first[_highest_active])
+ _highest_active--;
+ }
+ }
+
+ ///@}
+
+ ///Lift an active item to a higher level.
+
+ ///Lift an active item to a higher level.
+ ///\param i The item to be lifted. It must be active.
+ ///\param new_level The new level of \c i. It must be strictly higher
+ ///than the current level.
+ ///
+ void lift(Item i, int new_level)
+ {
+ const int lo = _level[i];
+ const Vit w = _where[i];
+
+ copy(_last_active[lo],w);
+ copy(--_first[lo+1],_last_active[lo]--);
+ for(int l=lo+1;l_highest_active) _highest_active=new_level;
+ }
+
+ ///Move an inactive item to the top but one level (in a dirty way).
+
+ ///This function moves an inactive item from the top level to the top
+ ///but one level (in a dirty way).
+ ///\warning It makes the underlying datastructure corrupt, so use it
+ ///only if you really know what it is for.
+ ///\pre The item is on the top level.
+ void dirtyTopButOne(Item i) {
+ _level.set(i,_max_level - 1);
+ }
+
+ ///Lift all items on and above the given level to the top level.
+
+ ///This function lifts all items on and above level \c l to the top
+ ///level and deactivates them.
+ void liftToTop(int l)
+ {
+ const Vit f=_first[l];
+ const Vit tl=_first[_max_level];
+ for(Vit i=f;i!=tl;++i)
+ _level.set(*i,_max_level);
+ for(int i=l;i<=_max_level;i++)
+ {
+ _first[i]=f;
+ _last_active[i]=f-1;
+ }
+ for(_highest_active=l-1;
+ _highest_active>=0 &&
+ _last_active[_highest_active]<_first[_highest_active];
+ _highest_active--) ;
+ }
+
+ private:
+ int _init_lev;
+ Vit _init_num;
+
+ public:
+
+ ///\name Initialization
+ ///Using these functions you can initialize the levels of the items.
+ ///\n
+ ///The initialization must be started with calling \c initStart().
+ ///Then the items should be listed level by level starting with the
+ ///lowest one (level 0) using \c initAddItem() and \c initNewLevel().
+ ///Finally \c initFinish() must be called.
+ ///The items not listed are put on the highest level.
+ ///@{
+
+ ///Start the initialization process.
+ void initStart()
+ {
+ _init_lev=0;
+ _init_num=&_items[0];
+ _first[0]=&_items[0];
+ _last_active[0]=&_items[0]-1;
+ Vit n=&_items[0];
+ for(typename ItemSetTraits::ItemIt i(_g);i!=INVALID;++i)
+ {
+ *n=i;
+ _where.set(i,n);
+ _level.set(i,_max_level);
+ ++n;
+ }
+ }
+
+ ///Add an item to the current level.
+ void initAddItem(Item i)
+ {
+ swap(_where[i],_init_num);
+ _level.set(i,_init_lev);
+ ++_init_num;
+ }
+
+ ///Start a new level.
+
+ ///Start a new level.
+ ///It shouldn't be used before the items on level 0 are listed.
+ void initNewLevel()
+ {
+ _init_lev++;
+ _first[_init_lev]=_init_num;
+ _last_active[_init_lev]=_init_num-1;
+ }
+
+ ///Finalize the initialization process.
+ void initFinish()
+ {
+ for(_init_lev++;_init_lev<=_max_level;_init_lev++)
+ {
+ _first[_init_lev]=_init_num;
+ _last_active[_init_lev]=_init_num-1;
+ }
+ _first[_max_level+1]=&_items[0]+_item_num;
+ _last_active[_max_level+1]=&_items[0]+_item_num-1;
+ _highest_active = -1;
+ }
+
+ ///@}
+
+ };
+
+ ///Class for handling "labels" in push-relabel type algorithms.
+
+ ///A class for handling "labels" in push-relabel type algorithms.
+ ///
+ ///\ingroup auxdat
+ ///Using this class you can assign "labels" (nonnegative integer numbers)
+ ///to the edges or nodes of a graph, manipulate and query them through
+ ///operations typically arising in "push-relabel" type algorithms.
+ ///
+ ///Each item is either \em active or not, and you can also choose a
+ ///highest level active item.
+ ///
+ ///\sa Elevator
+ ///
+ ///\param Graph Type of the underlying graph.
+ ///\param Item Type of the items the data is assigned to (Graph::Node,
+ ///Graph::Arc, Graph::Edge).
+ template
+ class LinkedElevator {
+ public:
+
+ typedef Item Key;
+ typedef int Value;
+
+ private:
+
+ typedef typename ItemSetTraits::
+ template Map
- ::Type ItemMap;
+ typedef typename ItemSetTraits::
+ template Map::Type IntMap;
+ typedef typename ItemSetTraits::
+ template Map::Type BoolMap;
+
+ const Graph &_graph;
+ int _max_level;
+ int _item_num;
+ std::vector
- _first, _last;
+ ItemMap _prev, _next;
+ int _highest_active;
+ IntMap _level;
+ BoolMap _active;
+
+ public:
+ ///Constructor with given maximum level.
+
+ ///Constructor with given maximum level.
+ ///
+ ///\param graph The underlying graph.
+ ///\param max_level The maximum allowed level.
+ ///Set the range of the possible labels to [0..max_level].
+ LinkedElevator(const Graph& graph, int max_level)
+ : _graph(graph), _max_level(max_level), _item_num(_max_level),
+ _first(_max_level + 1), _last(_max_level + 1),
+ _prev(graph), _next(graph),
+ _highest_active(-1), _level(graph), _active(graph) {}
+
+ ///Constructor.
+
+ ///Constructor.
+ ///
+ ///\param graph The underlying graph.
+ ///Set the range of the possible labels to [0..max_level],
+ ///where \c max_level is equal to the number of labeled items in the graph.
+ LinkedElevator(const Graph& graph)
+ : _graph(graph), _max_level(countItems(graph)),
+ _item_num(_max_level),
+ _first(_max_level + 1), _last(_max_level + 1),
+ _prev(graph, INVALID), _next(graph, INVALID),
+ _highest_active(-1), _level(graph), _active(graph) {}
+
+
+ ///Activate item \c i.
+
+ ///Activate item \c i.
+ ///\pre Item \c i shouldn't be active before.
+ void activate(Item i) {
+ _active.set(i, true);
+
+ int level = _level[i];
+ if (level > _highest_active) {
+ _highest_active = level;
+ }
+
+ if (_prev[i] == INVALID || _active[_prev[i]]) return;
+ //unlace
+ _next.set(_prev[i], _next[i]);
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], _prev[i]);
+ } else {
+ _last[level] = _prev[i];
+ }
+ //lace
+ _next.set(i, _first[level]);
+ _prev.set(_first[level], i);
+ _prev.set(i, INVALID);
+ _first[level] = i;
+
+ }
+
+ ///Deactivate item \c i.
+
+ ///Deactivate item \c i.
+ ///\pre Item \c i must be active before.
+ void deactivate(Item i) {
+ _active.set(i, false);
+ int level = _level[i];
+
+ if (_next[i] == INVALID || !_active[_next[i]])
+ goto find_highest_level;
+
+ //unlace
+ _prev.set(_next[i], _prev[i]);
+ if (_prev[i] != INVALID) {
+ _next.set(_prev[i], _next[i]);
+ } else {
+ _first[_level[i]] = _next[i];
+ }
+ //lace
+ _prev.set(i, _last[level]);
+ _next.set(_last[level], i);
+ _next.set(i, INVALID);
+ _last[level] = i;
+
+ find_highest_level:
+ if (level == _highest_active) {
+ while (_highest_active >= 0 && activeFree(_highest_active))
+ --_highest_active;
+ }
+ }
+
+ ///Query whether item \c i is active
+ bool active(Item i) const { return _active[i]; }
+
+ ///Return the level of item \c i.
+ int operator[](Item i) const { return _level[i]; }
+
+ ///Return the number of items on level \c l.
+ int onLevel(int l) const {
+ int num = 0;
+ Item n = _first[l];
+ while (n != INVALID) {
+ ++num;
+ n = _next[n];
+ }
+ return num;
+ }
+
+ ///Return true if the level is empty.
+ bool emptyLevel(int l) const {
+ return _first[l] == INVALID;
+ }
+
+ ///Return the number of items above level \c l.
+ int aboveLevel(int l) const {
+ int num = 0;
+ for (int level = l + 1; level < _max_level; ++level)
+ num += onLevel(level);
+ return num;
+ }
+
+ ///Return the number of active items on level \c l.
+ int activesOnLevel(int l) const {
+ int num = 0;
+ Item n = _first[l];
+ while (n != INVALID && _active[n]) {
+ ++num;
+ n = _next[n];
+ }
+ return num;
+ }
+
+ ///Return true if there is no active item on level \c l.
+ bool activeFree(int l) const {
+ return _first[l] == INVALID || !_active[_first[l]];
+ }
+
+ ///Return the maximum allowed level.
+ int maxLevel() const {
+ return _max_level;
+ }
+
+ ///\name Highest Active Item
+ ///Functions for working with the highest level
+ ///active item.
+
+ ///@{
+
+ ///Return a highest level active item.
+
+ ///Return a highest level active item or INVALID if there is no active
+ ///item.
+ Item highestActive() const {
+ return _highest_active >= 0 ? _first[_highest_active] : INVALID;
+ }
+
+ ///Return the highest active level.
+
+ ///Return the level of the highest active item or -1 if there is no active
+ ///item.
+ int highestActiveLevel() const {
+ return _highest_active;
+ }
+
+ ///Lift the highest active item by one.
+
+ ///Lift the item returned by highestActive() by one.
+ ///
+ void liftHighestActive() {
+ Item i = _first[_highest_active];
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], INVALID);
+ _first[_highest_active] = _next[i];
+ } else {
+ _first[_highest_active] = INVALID;
+ _last[_highest_active] = INVALID;
+ }
+ _level.set(i, ++_highest_active);
+ if (_first[_highest_active] == INVALID) {
+ _first[_highest_active] = i;
+ _last[_highest_active] = i;
+ _prev.set(i, INVALID);
+ _next.set(i, INVALID);
+ } else {
+ _prev.set(_first[_highest_active], i);
+ _next.set(i, _first[_highest_active]);
+ _first[_highest_active] = i;
+ }
+ }
+
+ ///Lift the highest active item to the given level.
+
+ ///Lift the item returned by highestActive() to level \c new_level.
+ ///
+ ///\warning \c new_level must be strictly higher
+ ///than the current level.
+ ///
+ void liftHighestActive(int new_level) {
+ Item i = _first[_highest_active];
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], INVALID);
+ _first[_highest_active] = _next[i];
+ } else {
+ _first[_highest_active] = INVALID;
+ _last[_highest_active] = INVALID;
+ }
+ _level.set(i, _highest_active = new_level);
+ if (_first[_highest_active] == INVALID) {
+ _first[_highest_active] = _last[_highest_active] = i;
+ _prev.set(i, INVALID);
+ _next.set(i, INVALID);
+ } else {
+ _prev.set(_first[_highest_active], i);
+ _next.set(i, _first[_highest_active]);
+ _first[_highest_active] = i;
+ }
+ }
+
+ ///Lift the highest active item to the top level.
+
+ ///Lift the item returned by highestActive() to the top level and
+ ///deactivate it.
+ void liftHighestActiveToTop() {
+ Item i = _first[_highest_active];
+ _level.set(i, _max_level);
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], INVALID);
+ _first[_highest_active] = _next[i];
+ } else {
+ _first[_highest_active] = INVALID;
+ _last[_highest_active] = INVALID;
+ }
+ while (_highest_active >= 0 && activeFree(_highest_active))
+ --_highest_active;
+ }
+
+ ///@}
+
+ ///\name Active Item on Certain Level
+ ///Functions for working with the active items.
+
+ ///@{
+
+ ///Return an active item on level \c l.
+
+ ///Return an active item on level \c l or \ref INVALID if there is no such
+ ///an item. (\c l must be from the range [0...\c max_level].
+ Item activeOn(int l) const
+ {
+ return _active[_first[l]] ? _first[l] : INVALID;
+ }
+
+ ///Lift the active item returned by \c activeOn(l) by one.
+
+ ///Lift the active item returned by \ref activeOn() "activeOn(l)"
+ ///by one.
+ Item liftActiveOn(int l)
+ {
+ Item i = _first[l];
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], INVALID);
+ _first[l] = _next[i];
+ } else {
+ _first[l] = INVALID;
+ _last[l] = INVALID;
+ }
+ _level.set(i, ++l);
+ if (_first[l] == INVALID) {
+ _first[l] = _last[l] = i;
+ _prev.set(i, INVALID);
+ _next.set(i, INVALID);
+ } else {
+ _prev.set(_first[l], i);
+ _next.set(i, _first[l]);
+ _first[l] = i;
+ }
+ if (_highest_active < l) {
+ _highest_active = l;
+ }
+ }
+
+ ///Lift the active item returned by \c activeOn(l) to the given level.
+
+ ///Lift the active item returned by \ref activeOn() "activeOn(l)"
+ ///to the given level.
+ void liftActiveOn(int l, int new_level)
+ {
+ Item i = _first[l];
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], INVALID);
+ _first[l] = _next[i];
+ } else {
+ _first[l] = INVALID;
+ _last[l] = INVALID;
+ }
+ _level.set(i, l = new_level);
+ if (_first[l] == INVALID) {
+ _first[l] = _last[l] = i;
+ _prev.set(i, INVALID);
+ _next.set(i, INVALID);
+ } else {
+ _prev.set(_first[l], i);
+ _next.set(i, _first[l]);
+ _first[l] = i;
+ }
+ if (_highest_active < l) {
+ _highest_active = l;
+ }
+ }
+
+ ///Lift the active item returned by \c activeOn(l) to the top level.
+
+ ///Lift the active item returned by \ref activeOn() "activeOn(l)"
+ ///to the top level and deactivate it.
+ void liftActiveToTop(int l)
+ {
+ Item i = _first[l];
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], INVALID);
+ _first[l] = _next[i];
+ } else {
+ _first[l] = INVALID;
+ _last[l] = INVALID;
+ }
+ _level.set(i, _max_level);
+ if (l == _highest_active) {
+ while (_highest_active >= 0 && activeFree(_highest_active))
+ --_highest_active;
+ }
+ }
+
+ ///@}
+
+ /// \brief Lift an active item to a higher level.
+ ///
+ /// Lift an active item to a higher level.
+ /// \param i The item to be lifted. It must be active.
+ /// \param new_level The new level of \c i. It must be strictly higher
+ /// than the current level.
+ ///
+ void lift(Item i, int new_level) {
+ if (_next[i] != INVALID) {
+ _prev.set(_next[i], _prev[i]);
+ } else {
+ _last[new_level] = _prev[i];
+ }
+ if (_prev[i] != INVALID) {
+ _next.set(_prev[i], _next[i]);
+ } else {
+ _first[new_level] = _next[i];
+ }
+ _level.set(i, new_level);
+ if (_first[new_level] == INVALID) {
+ _first[new_level] = _last[new_level] = i;
+ _prev.set(i, INVALID);
+ _next.set(i, INVALID);
+ } else {
+ _prev.set(_first[new_level], i);
+ _next.set(i, _first[new_level]);
+ _first[new_level] = i;
+ }
+ if (_highest_active < new_level) {
+ _highest_active = new_level;
+ }
+ }
+
+ ///Move an inactive item to the top but one level (in a dirty way).
+
+ ///This function moves an inactive item from the top level to the top
+ ///but one level (in a dirty way).
+ ///\warning It makes the underlying datastructure corrupt, so use it
+ ///only if you really know what it is for.
+ ///\pre The item is on the top level.
+ void dirtyTopButOne(Item i) {
+ _level.set(i, _max_level - 1);
+ }
+
+ ///Lift all items on and above the given level to the top level.
+
+ ///This function lifts all items on and above level \c l to the top
+ ///level and deactivates them.
+ void liftToTop(int l) {
+ for (int i = l + 1; _first[i] != INVALID; ++i) {
+ Item n = _first[i];
+ while (n != INVALID) {
+ _level.set(n, _max_level);
+ n = _next[n];
+ }
+ _first[i] = INVALID;
+ _last[i] = INVALID;
+ }
+ if (_highest_active > l - 1) {
+ _highest_active = l - 1;
+ while (_highest_active >= 0 && activeFree(_highest_active))
+ --_highest_active;
+ }
+ }
+
+ private:
+
+ int _init_level;
+
+ public:
+
+ ///\name Initialization
+ ///Using these functions you can initialize the levels of the items.
+ ///\n
+ ///The initialization must be started with calling \c initStart().
+ ///Then the items should be listed level by level starting with the
+ ///lowest one (level 0) using \c initAddItem() and \c initNewLevel().
+ ///Finally \c initFinish() must be called.
+ ///The items not listed are put on the highest level.
+ ///@{
+
+ ///Start the initialization process.
+ void initStart() {
+
+ for (int i = 0; i <= _max_level; ++i) {
+ _first[i] = _last[i] = INVALID;
+ }
+ _init_level = 0;
+ for(typename ItemSetTraits::ItemIt i(_graph);
+ i != INVALID; ++i) {
+ _level.set(i, _max_level);
+ _active.set(i, false);
+ }
+ }
+
+ ///Add an item to the current level.
+ void initAddItem(Item i) {
+ _level.set(i, _init_level);
+ if (_last[_init_level] == INVALID) {
+ _first[_init_level] = i;
+ _last[_init_level] = i;
+ _prev.set(i, INVALID);
+ _next.set(i, INVALID);
+ } else {
+ _prev.set(i, _last[_init_level]);
+ _next.set(i, INVALID);
+ _next.set(_last[_init_level], i);
+ _last[_init_level] = i;
+ }
+ }
+
+ ///Start a new level.
+
+ ///Start a new level.
+ ///It shouldn't be used before the items on level 0 are listed.
+ void initNewLevel() {
+ ++_init_level;
+ }
+
+ ///Finalize the initialization process.
+ void initFinish() {
+ _highest_active = -1;
+ }
+
+ ///@}
+
+ };
+
+
+} //END OF NAMESPACE LEMON
+
+#endif
+
Index: lemon/full_graph.h
===================================================================
--- lemon/full_graph.h (revision 372)
+++ lemon/full_graph.h (revision 372)
@@ -0,0 +1,614 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_FULL_GRAPH_H
+#define LEMON_FULL_GRAPH_H
+
+#include
+#include
+
+///\ingroup graphs
+///\file
+///\brief FullGraph and FullDigraph classes.
+
+namespace lemon {
+
+ class FullDigraphBase {
+ public:
+
+ typedef FullDigraphBase Graph;
+
+ class Node;
+ class Arc;
+
+ protected:
+
+ int _node_num;
+ int _arc_num;
+
+ FullDigraphBase() {}
+
+ void construct(int n) { _node_num = n; _arc_num = n * n; }
+
+ public:
+
+ typedef True NodeNumTag;
+ typedef True ArcNumTag;
+
+ Node operator()(int ix) const { return Node(ix); }
+ int index(const Node& node) const { return node._id; }
+
+ Arc arc(const Node& s, const Node& t) const {
+ return Arc(s._id * _node_num + t._id);
+ }
+
+ int nodeNum() const { return _node_num; }
+ int arcNum() const { return _arc_num; }
+
+ int maxNodeId() const { return _node_num - 1; }
+ int maxArcId() const { return _arc_num - 1; }
+
+ Node source(Arc arc) const { return arc._id / _node_num; }
+ Node target(Arc arc) const { return arc._id % _node_num; }
+
+ static int id(Node node) { return node._id; }
+ static int id(Arc arc) { return arc._id; }
+
+ static Node nodeFromId(int id) { return Node(id);}
+ static Arc arcFromId(int id) { return Arc(id);}
+
+ typedef True FindArcTag;
+
+ Arc findArc(Node s, Node t, Arc prev = INVALID) const {
+ return prev == INVALID ? arc(s, t) : INVALID;
+ }
+
+ class Node {
+ friend class FullDigraphBase;
+
+ protected:
+ int _id;
+ Node(int id) : _id(id) {}
+ public:
+ Node() {}
+ Node (Invalid) : _id(-1) {}
+ bool operator==(const Node node) const {return _id == node._id;}
+ bool operator!=(const Node node) const {return _id != node._id;}
+ bool operator<(const Node node) const {return _id < node._id;}
+ };
+
+ class Arc {
+ friend class FullDigraphBase;
+
+ protected:
+ int _id; // _node_num * source + target;
+
+ Arc(int id) : _id(id) {}
+
+ public:
+ Arc() { }
+ Arc (Invalid) { _id = -1; }
+ bool operator==(const Arc arc) const {return _id == arc._id;}
+ bool operator!=(const Arc arc) const {return _id != arc._id;}
+ bool operator<(const Arc arc) const {return _id < arc._id;}
+ };
+
+ void first(Node& node) const {
+ node._id = _node_num - 1;
+ }
+
+ static void next(Node& node) {
+ --node._id;
+ }
+
+ void first(Arc& arc) const {
+ arc._id = _arc_num - 1;
+ }
+
+ static void next(Arc& arc) {
+ --arc._id;
+ }
+
+ void firstOut(Arc& arc, const Node& node) const {
+ arc._id = (node._id + 1) * _node_num - 1;
+ }
+
+ void nextOut(Arc& arc) const {
+ if (arc._id % _node_num == 0) arc._id = 0;
+ --arc._id;
+ }
+
+ void firstIn(Arc& arc, const Node& node) const {
+ arc._id = _arc_num + node._id - _node_num;
+ }
+
+ void nextIn(Arc& arc) const {
+ arc._id -= _node_num;
+ if (arc._id < 0) arc._id = -1;
+ }
+
+ };
+
+ typedef DigraphExtender ExtendedFullDigraphBase;
+
+ /// \ingroup graphs
+ ///
+ /// \brief A full digraph class.
+ ///
+ /// This is a simple and fast directed full graph implementation.
+ /// From each node go arcs to each node (including the source node),
+ /// therefore the number of the arcs in the digraph is the square of
+ /// the node number. This digraph type is completely static, so you
+ /// can neither add nor delete either arcs or nodes, and it needs
+ /// constant space in memory.
+ ///
+ /// This class conforms to the \ref concepts::Digraph "Digraph" concept
+ /// and it also has an important extra feature that its maps are
+ /// real \ref concepts::ReferenceMap "reference map"s.
+ ///
+ /// The \c FullDigraph and \c FullGraph classes are very similar,
+ /// but there are two differences. While this class conforms only
+ /// to the \ref concepts::Digraph "Digraph" concept, the \c FullGraph
+ /// class conforms to the \ref concepts::Graph "Graph" concept,
+ /// moreover \c FullGraph does not contain a loop arc for each
+ /// node as \c FullDigraph does.
+ ///
+ /// \sa FullGraph
+ class FullDigraph : public ExtendedFullDigraphBase {
+ public:
+
+ typedef ExtendedFullDigraphBase Parent;
+
+ /// \brief Constructor
+ FullDigraph() { construct(0); }
+
+ /// \brief Constructor
+ ///
+ /// Constructor.
+ /// \param n The number of the nodes.
+ FullDigraph(int n) { construct(n); }
+
+ /// \brief Resizes the digraph
+ ///
+ /// Resizes the digraph. The function will fully destroy and
+ /// rebuild the digraph. This cause that the maps of the digraph will
+ /// reallocated automatically and the previous values will be lost.
+ void resize(int n) {
+ Parent::notifier(Arc()).clear();
+ Parent::notifier(Node()).clear();
+ construct(n);
+ Parent::notifier(Node()).build();
+ Parent::notifier(Arc()).build();
+ }
+
+ /// \brief Returns the node with the given index.
+ ///
+ /// Returns the node with the given index. Since it is a static
+ /// digraph its nodes can be indexed with integers from the range
+ /// [0..nodeNum()-1].
+ /// \sa index()
+ Node operator()(int ix) const { return Parent::operator()(ix); }
+
+ /// \brief Returns the index of the given node.
+ ///
+ /// Returns the index of the given node. Since it is a static
+ /// digraph its nodes can be indexed with integers from the range
+ /// [0..nodeNum()-1].
+ /// \sa operator()
+ int index(const Node& node) const { return Parent::index(node); }
+
+ /// \brief Returns the arc connecting the given nodes.
+ ///
+ /// Returns the arc connecting the given nodes.
+ Arc arc(const Node& u, const Node& v) const {
+ return Parent::arc(u, v);
+ }
+
+ /// \brief Number of nodes.
+ int nodeNum() const { return Parent::nodeNum(); }
+ /// \brief Number of arcs.
+ int arcNum() const { return Parent::arcNum(); }
+ };
+
+
+ class FullGraphBase {
+ int _node_num;
+ int _edge_num;
+ public:
+
+ typedef FullGraphBase Graph;
+
+ class Node;
+ class Arc;
+ class Edge;
+
+ protected:
+
+ FullGraphBase() {}
+
+ void construct(int n) { _node_num = n; _edge_num = n * (n - 1) / 2; }
+
+ int _uid(int e) const {
+ int u = e / _node_num;
+ int v = e % _node_num;
+ return u < v ? u : _node_num - 2 - u;
+ }
+
+ int _vid(int e) const {
+ int u = e / _node_num;
+ int v = e % _node_num;
+ return u < v ? v : _node_num - 1 - v;
+ }
+
+ void _uvid(int e, int& u, int& v) const {
+ u = e / _node_num;
+ v = e % _node_num;
+ if (u >= v) {
+ u = _node_num - 2 - u;
+ v = _node_num - 1 - v;
+ }
+ }
+
+ void _stid(int a, int& s, int& t) const {
+ if ((a & 1) == 1) {
+ _uvid(a >> 1, s, t);
+ } else {
+ _uvid(a >> 1, t, s);
+ }
+ }
+
+ int _eid(int u, int v) const {
+ if (u < (_node_num - 1) / 2) {
+ return u * _node_num + v;
+ } else {
+ return (_node_num - 1 - u) * _node_num - v - 1;
+ }
+ }
+
+ public:
+
+ Node operator()(int ix) const { return Node(ix); }
+ int index(const Node& node) const { return node._id; }
+
+ Edge edge(const Node& u, const Node& v) const {
+ if (u._id < v._id) {
+ return Edge(_eid(u._id, v._id));
+ } else if (u._id != v._id) {
+ return Edge(_eid(v._id, u._id));
+ } else {
+ return INVALID;
+ }
+ }
+
+ Arc arc(const Node& s, const Node& t) const {
+ if (s._id < t._id) {
+ return Arc((_eid(s._id, t._id) << 1) | 1);
+ } else if (s._id != t._id) {
+ return Arc(_eid(t._id, s._id) << 1);
+ } else {
+ return INVALID;
+ }
+ }
+
+ typedef True NodeNumTag;
+ typedef True ArcNumTag;
+ typedef True EdgeNumTag;
+
+ int nodeNum() const { return _node_num; }
+ int arcNum() const { return 2 * _edge_num; }
+ int edgeNum() const { return _edge_num; }
+
+ static int id(Node node) { return node._id; }
+ static int id(Arc arc) { return arc._id; }
+ static int id(Edge edge) { return edge._id; }
+
+ int maxNodeId() const { return _node_num-1; }
+ int maxArcId() const { return 2 * _edge_num-1; }
+ int maxEdgeId() const { return _edge_num-1; }
+
+ static Node nodeFromId(int id) { return Node(id);}
+ static Arc arcFromId(int id) { return Arc(id);}
+ static Edge edgeFromId(int id) { return Edge(id);}
+
+ Node u(Edge edge) const {
+ return Node(_uid(edge._id));
+ }
+
+ Node v(Edge edge) const {
+ return Node(_vid(edge._id));
+ }
+
+ Node source(Arc arc) const {
+ return Node((arc._id & 1) == 1 ?
+ _uid(arc._id >> 1) : _vid(arc._id >> 1));
+ }
+
+ Node target(Arc arc) const {
+ return Node((arc._id & 1) == 1 ?
+ _vid(arc._id >> 1) : _uid(arc._id >> 1));
+ }
+
+ typedef True FindEdgeTag;
+ typedef True FindArcTag;
+
+ Edge findEdge(Node u, Node v, Edge prev = INVALID) const {
+ return prev != INVALID ? INVALID : edge(u, v);
+ }
+
+ Arc findArc(Node s, Node t, Arc prev = INVALID) const {
+ return prev != INVALID ? INVALID : arc(s, t);
+ }
+
+ class Node {
+ friend class FullGraphBase;
+
+ protected:
+ int _id;
+ Node(int id) : _id(id) {}
+ public:
+ Node() {}
+ Node (Invalid) { _id = -1; }
+ bool operator==(const Node node) const {return _id == node._id;}
+ bool operator!=(const Node node) const {return _id != node._id;}
+ bool operator<(const Node node) const {return _id < node._id;}
+ };
+
+ class Edge {
+ friend class FullGraphBase;
+ friend class Arc;
+
+ protected:
+ int _id;
+
+ Edge(int id) : _id(id) {}
+
+ public:
+ Edge() { }
+ Edge (Invalid) { _id = -1; }
+
+ bool operator==(const Edge edge) const {return _id == edge._id;}
+ bool operator!=(const Edge edge) const {return _id != edge._id;}
+ bool operator<(const Edge edge) const {return _id < edge._id;}
+ };
+
+ class Arc {
+ friend class FullGraphBase;
+
+ protected:
+ int _id;
+
+ Arc(int id) : _id(id) {}
+
+ public:
+ Arc() { }
+ Arc (Invalid) { _id = -1; }
+
+ operator Edge() const { return Edge(_id != -1 ? (_id >> 1) : -1); }
+
+ bool operator==(const Arc arc) const {return _id == arc._id;}
+ bool operator!=(const Arc arc) const {return _id != arc._id;}
+ bool operator<(const Arc arc) const {return _id < arc._id;}
+ };
+
+ static bool direction(Arc arc) {
+ return (arc._id & 1) == 1;
+ }
+
+ static Arc direct(Edge edge, bool dir) {
+ return Arc((edge._id << 1) | (dir ? 1 : 0));
+ }
+
+ void first(Node& node) const {
+ node._id = _node_num - 1;
+ }
+
+ static void next(Node& node) {
+ --node._id;
+ }
+
+ void first(Arc& arc) const {
+ arc._id = (_edge_num << 1) - 1;
+ }
+
+ static void next(Arc& arc) {
+ --arc._id;
+ }
+
+ void first(Edge& edge) const {
+ edge._id = _edge_num - 1;
+ }
+
+ static void next(Edge& edge) {
+ --edge._id;
+ }
+
+ void firstOut(Arc& arc, const Node& node) const {
+ int s = node._id, t = _node_num - 1;
+ if (s < t) {
+ arc._id = (_eid(s, t) << 1) | 1;
+ } else {
+ --t;
+ arc._id = (t != -1 ? (_eid(t, s) << 1) : -1);
+ }
+ }
+
+ void nextOut(Arc& arc) const {
+ int s, t;
+ _stid(arc._id, s, t);
+ --t;
+ if (s < t) {
+ arc._id = (_eid(s, t) << 1) | 1;
+ } else {
+ if (s == t) --t;
+ arc._id = (t != -1 ? (_eid(t, s) << 1) : -1);
+ }
+ }
+
+ void firstIn(Arc& arc, const Node& node) const {
+ int s = _node_num - 1, t = node._id;
+ if (s > t) {
+ arc._id = (_eid(t, s) << 1);
+ } else {
+ --s;
+ arc._id = (s != -1 ? (_eid(s, t) << 1) | 1 : -1);
+ }
+ }
+
+ void nextIn(Arc& arc) const {
+ int s, t;
+ _stid(arc._id, s, t);
+ --s;
+ if (s > t) {
+ arc._id = (_eid(t, s) << 1);
+ } else {
+ if (s == t) --s;
+ arc._id = (s != -1 ? (_eid(s, t) << 1) | 1 : -1);
+ }
+ }
+
+ void firstInc(Edge& edge, bool& dir, const Node& node) const {
+ int u = node._id, v = _node_num - 1;
+ if (u < v) {
+ edge._id = _eid(u, v);
+ dir = true;
+ } else {
+ --v;
+ edge._id = (v != -1 ? _eid(v, u) : -1);
+ dir = false;
+ }
+ }
+
+ void nextInc(Edge& edge, bool& dir) const {
+ int u, v;
+ if (dir) {
+ _uvid(edge._id, u, v);
+ --v;
+ if (u < v) {
+ edge._id = _eid(u, v);
+ } else {
+ --v;
+ edge._id = (v != -1 ? _eid(v, u) : -1);
+ dir = false;
+ }
+ } else {
+ _uvid(edge._id, v, u);
+ --v;
+ edge._id = (v != -1 ? _eid(v, u) : -1);
+ }
+ }
+
+ };
+
+ typedef GraphExtender ExtendedFullGraphBase;
+
+ /// \ingroup graphs
+ ///
+ /// \brief An undirected full graph class.
+ ///
+ /// This is a simple and fast undirected full graph
+ /// implementation. From each node go edge to each other node,
+ /// therefore the number of edges in the graph is \f$n(n-1)/2\f$.
+ /// This graph type is completely static, so you can neither
+ /// add nor delete either edges or nodes, and it needs constant
+ /// space in memory.
+ ///
+ /// This class conforms to the \ref concepts::Graph "Graph" concept
+ /// and it also has an important extra feature that its maps are
+ /// real \ref concepts::ReferenceMap "reference map"s.
+ ///
+ /// The \c FullGraph and \c FullDigraph classes are very similar,
+ /// but there are two differences. While the \c FullDigraph class
+ /// conforms only to the \ref concepts::Digraph "Digraph" concept,
+ /// this class conforms to the \ref concepts::Graph "Graph" concept,
+ /// moreover \c FullGraph does not contain a loop arc for each
+ /// node as \c FullDigraph does.
+ ///
+ /// \sa FullDigraph
+ class FullGraph : public ExtendedFullGraphBase {
+ public:
+
+ typedef ExtendedFullGraphBase Parent;
+
+ /// \brief Constructor
+ FullGraph() { construct(0); }
+
+ /// \brief Constructor
+ ///
+ /// Constructor.
+ /// \param n The number of the nodes.
+ FullGraph(int n) { construct(n); }
+
+ /// \brief Resizes the graph
+ ///
+ /// Resizes the graph. The function will fully destroy and
+ /// rebuild the graph. This cause that the maps of the graph will
+ /// reallocated automatically and the previous values will be lost.
+ void resize(int n) {
+ Parent::notifier(Arc()).clear();
+ Parent::notifier(Edge()).clear();
+ Parent::notifier(Node()).clear();
+ construct(n);
+ Parent::notifier(Node()).build();
+ Parent::notifier(Edge()).build();
+ Parent::notifier(Arc()).build();
+ }
+
+ /// \brief Returns the node with the given index.
+ ///
+ /// Returns the node with the given index. Since it is a static
+ /// graph its nodes can be indexed with integers from the range
+ /// [0..nodeNum()-1].
+ /// \sa index()
+ Node operator()(int ix) const { return Parent::operator()(ix); }
+
+ /// \brief Returns the index of the given node.
+ ///
+ /// Returns the index of the given node. Since it is a static
+ /// graph its nodes can be indexed with integers from the range
+ /// [0..nodeNum()-1].
+ /// \sa operator()
+ int index(const Node& node) const { return Parent::index(node); }
+
+ /// \brief Returns the arc connecting the given nodes.
+ ///
+ /// Returns the arc connecting the given nodes.
+ Arc arc(const Node& s, const Node& t) const {
+ return Parent::arc(s, t);
+ }
+
+ /// \brief Returns the edge connects the given nodes.
+ ///
+ /// Returns the edge connects the given nodes.
+ Edge edge(const Node& u, const Node& v) const {
+ return Parent::edge(u, v);
+ }
+
+ /// \brief Number of nodes.
+ int nodeNum() const { return Parent::nodeNum(); }
+ /// \brief Number of arcs.
+ int arcNum() const { return Parent::arcNum(); }
+ /// \brief Number of edges.
+ int edgeNum() const { return Parent::edgeNum(); }
+
+ };
+
+
+} //namespace lemon
+
+
+#endif //LEMON_FULL_GRAPH_H
Index: lemon/grid_graph.h
===================================================================
--- lemon/grid_graph.h (revision 372)
+++ lemon/grid_graph.h (revision 372)
@@ -0,0 +1,705 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef GRID_GRAPH_H
+#define GRID_GRAPH_H
+
+#include
+#include
+#include
+#include
+
+///\ingroup graphs
+///\file
+///\brief GridGraph class.
+
+namespace lemon {
+
+ class GridGraphBase {
+
+ public:
+
+ typedef GridGraphBase Graph;
+
+ class Node;
+ class Edge;
+ class Arc;
+
+ public:
+
+ GridGraphBase() {}
+
+ protected:
+
+ void construct(int width, int height) {
+ _width = width; _height = height;
+ _node_num = width * height;
+ _edge_num = 2 * _node_num - width - height;
+ _edge_limit = _node_num - _width;
+ }
+
+ public:
+
+ Node operator()(int i, int j) const {
+ LEMON_DEBUG(0 <= i && i < _width &&
+ 0 <= j && j < _height, "Index out of range");
+ return Node(i + j * _width);
+ }
+
+ int col(Node n) const {
+ return n._id % _width;
+ }
+
+ int row(Node n) const {
+ return n._id / _width;
+ }
+
+ dim2::Point pos(Node n) const {
+ return dim2::Point(col(n), row(n));
+ }
+
+ int width() const {
+ return _width;
+ }
+
+ int height() const {
+ return _height;
+ }
+
+ typedef True NodeNumTag;
+ typedef True EdgeNumTag;
+ typedef True ArcNumTag;
+
+ int nodeNum() const { return _node_num; }
+ int edgeNum() const { return _edge_num; }
+ int arcNum() const { return 2 * _edge_num; }
+
+ Node u(Edge edge) const {
+ if (edge._id < _edge_limit) {
+ return edge._id;
+ } else {
+ return (edge._id - _edge_limit) % (_width - 1) +
+ (edge._id - _edge_limit) / (_width - 1) * _width;
+ }
+ }
+
+ Node v(Edge edge) const {
+ if (edge._id < _edge_limit) {
+ return edge._id + _width;
+ } else {
+ return (edge._id - _edge_limit) % (_width - 1) +
+ (edge._id - _edge_limit) / (_width - 1) * _width + 1;
+ }
+ }
+
+ Node source(Arc arc) const {
+ return (arc._id & 1) == 1 ? u(arc) : v(arc);
+ }
+
+ Node target(Arc arc) const {
+ return (arc._id & 1) == 1 ? v(arc) : u(arc);
+ }
+
+ static int id(Node node) { return node._id; }
+ static int id(Edge edge) { return edge._id; }
+ static int id(Arc arc) { return arc._id; }
+
+ int maxNodeId() const { return _node_num - 1; }
+ int maxEdgeId() const { return _edge_num - 1; }
+ int maxArcId() const { return 2 * _edge_num - 1; }
+
+ static Node nodeFromId(int id) { return Node(id);}
+ static Edge edgeFromId(int id) { return Edge(id);}
+ static Arc arcFromId(int id) { return Arc(id);}
+
+ typedef True FindEdgeTag;
+ typedef True FindArcTag;
+
+ Edge findEdge(Node u, Node v, Edge prev = INVALID) const {
+ if (prev != INVALID) return INVALID;
+ if (v._id > u._id) {
+ if (v._id - u._id == _width)
+ return Edge(u._id);
+ if (v._id - u._id == 1 && u._id % _width < _width - 1) {
+ return Edge(u._id / _width * (_width - 1) +
+ u._id % _width + _edge_limit);
+ }
+ } else {
+ if (u._id - v._id == _width)
+ return Edge(v._id);
+ if (u._id - v._id == 1 && v._id % _width < _width - 1) {
+ return Edge(v._id / _width * (_width - 1) +
+ v._id % _width + _edge_limit);
+ }
+ }
+ return INVALID;
+ }
+
+ Arc findArc(Node u, Node v, Arc prev = INVALID) const {
+ if (prev != INVALID) return INVALID;
+ if (v._id > u._id) {
+ if (v._id - u._id == _width)
+ return Arc((u._id << 1) | 1);
+ if (v._id - u._id == 1 && u._id % _width < _width - 1) {
+ return Arc(((u._id / _width * (_width - 1) +
+ u._id % _width + _edge_limit) << 1) | 1);
+ }
+ } else {
+ if (u._id - v._id == _width)
+ return Arc(v._id << 1);
+ if (u._id - v._id == 1 && v._id % _width < _width - 1) {
+ return Arc((v._id / _width * (_width - 1) +
+ v._id % _width + _edge_limit) << 1);
+ }
+ }
+ return INVALID;
+ }
+
+ class Node {
+ friend class GridGraphBase;
+
+ protected:
+ int _id;
+ Node(int id) : _id(id) {}
+ public:
+ Node() {}
+ Node (Invalid) : _id(-1) {}
+ bool operator==(const Node node) const {return _id == node._id;}
+ bool operator!=(const Node node) const {return _id != node._id;}
+ bool operator<(const Node node) const {return _id < node._id;}
+ };
+
+ class Edge {
+ friend class GridGraphBase;
+ friend class Arc;
+
+ protected:
+ int _id;
+
+ Edge(int id) : _id(id) {}
+
+ public:
+ Edge() {}
+ Edge (Invalid) : _id(-1) {}
+ bool operator==(const Edge edge) const {return _id == edge._id;}
+ bool operator!=(const Edge edge) const {return _id != edge._id;}
+ bool operator<(const Edge edge) const {return _id < edge._id;}
+ };
+
+ class Arc {
+ friend class GridGraphBase;
+
+ protected:
+ int _id;
+
+ Arc(int id) : _id(id) {}
+
+ public:
+ Arc() {}
+ Arc (Invalid) : _id(-1) {}
+ operator Edge() const { return _id != -1 ? Edge(_id >> 1) : INVALID; }
+ bool operator==(const Arc arc) const {return _id == arc._id;}
+ bool operator!=(const Arc arc) const {return _id != arc._id;}
+ bool operator<(const Arc arc) const {return _id < arc._id;}
+ };
+
+ static bool direction(Arc arc) {
+ return (arc._id & 1) == 1;
+ }
+
+ static Arc direct(Edge edge, bool dir) {
+ return Arc((edge._id << 1) | (dir ? 1 : 0));
+ }
+
+ void first(Node& node) const {
+ node._id = _node_num - 1;
+ }
+
+ static void next(Node& node) {
+ --node._id;
+ }
+
+ void first(Edge& edge) const {
+ edge._id = _edge_num - 1;
+ }
+
+ static void next(Edge& edge) {
+ --edge._id;
+ }
+
+ void first(Arc& arc) const {
+ arc._id = 2 * _edge_num - 1;
+ }
+
+ static void next(Arc& arc) {
+ --arc._id;
+ }
+
+ void firstOut(Arc& arc, const Node& node) const {
+ if (node._id % _width < _width - 1) {
+ arc._id = (_edge_limit + node._id % _width +
+ (node._id / _width) * (_width - 1)) << 1 | 1;
+ return;
+ }
+ if (node._id < _node_num - _width) {
+ arc._id = node._id << 1 | 1;
+ return;
+ }
+ if (node._id % _width > 0) {
+ arc._id = (_edge_limit + node._id % _width +
+ (node._id / _width) * (_width - 1) - 1) << 1;
+ return;
+ }
+ if (node._id >= _width) {
+ arc._id = (node._id - _width) << 1;
+ return;
+ }
+ arc._id = -1;
+ }
+
+ void nextOut(Arc& arc) const {
+ int nid = arc._id >> 1;
+ if ((arc._id & 1) == 1) {
+ if (nid >= _edge_limit) {
+ nid = (nid - _edge_limit) % (_width - 1) +
+ (nid - _edge_limit) / (_width - 1) * _width;
+ if (nid < _node_num - _width) {
+ arc._id = nid << 1 | 1;
+ return;
+ }
+ }
+ if (nid % _width > 0) {
+ arc._id = (_edge_limit + nid % _width +
+ (nid / _width) * (_width - 1) - 1) << 1;
+ return;
+ }
+ if (nid >= _width) {
+ arc._id = (nid - _width) << 1;
+ return;
+ }
+ } else {
+ if (nid >= _edge_limit) {
+ nid = (nid - _edge_limit) % (_width - 1) +
+ (nid - _edge_limit) / (_width - 1) * _width + 1;
+ if (nid >= _width) {
+ arc._id = (nid - _width) << 1;
+ return;
+ }
+ }
+ }
+ arc._id = -1;
+ }
+
+ void firstIn(Arc& arc, const Node& node) const {
+ if (node._id % _width < _width - 1) {
+ arc._id = (_edge_limit + node._id % _width +
+ (node._id / _width) * (_width - 1)) << 1;
+ return;
+ }
+ if (node._id < _node_num - _width) {
+ arc._id = node._id << 1;
+ return;
+ }
+ if (node._id % _width > 0) {
+ arc._id = (_edge_limit + node._id % _width +
+ (node._id / _width) * (_width - 1) - 1) << 1 | 1;
+ return;
+ }
+ if (node._id >= _width) {
+ arc._id = (node._id - _width) << 1 | 1;
+ return;
+ }
+ arc._id = -1;
+ }
+
+ void nextIn(Arc& arc) const {
+ int nid = arc._id >> 1;
+ if ((arc._id & 1) == 0) {
+ if (nid >= _edge_limit) {
+ nid = (nid - _edge_limit) % (_width - 1) +
+ (nid - _edge_limit) / (_width - 1) * _width;
+ if (nid < _node_num - _width) {
+ arc._id = nid << 1;
+ return;
+ }
+ }
+ if (nid % _width > 0) {
+ arc._id = (_edge_limit + nid % _width +
+ (nid / _width) * (_width - 1) - 1) << 1 | 1;
+ return;
+ }
+ if (nid >= _width) {
+ arc._id = (nid - _width) << 1 | 1;
+ return;
+ }
+ } else {
+ if (nid >= _edge_limit) {
+ nid = (nid - _edge_limit) % (_width - 1) +
+ (nid - _edge_limit) / (_width - 1) * _width + 1;
+ if (nid >= _width) {
+ arc._id = (nid - _width) << 1 | 1;
+ return;
+ }
+ }
+ }
+ arc._id = -1;
+ }
+
+ void firstInc(Edge& edge, bool& dir, const Node& node) const {
+ if (node._id % _width < _width - 1) {
+ edge._id = _edge_limit + node._id % _width +
+ (node._id / _width) * (_width - 1);
+ dir = true;
+ return;
+ }
+ if (node._id < _node_num - _width) {
+ edge._id = node._id;
+ dir = true;
+ return;
+ }
+ if (node._id % _width > 0) {
+ edge._id = _edge_limit + node._id % _width +
+ (node._id / _width) * (_width - 1) - 1;
+ dir = false;
+ return;
+ }
+ if (node._id >= _width) {
+ edge._id = node._id - _width;
+ dir = false;
+ return;
+ }
+ edge._id = -1;
+ dir = true;
+ }
+
+ void nextInc(Edge& edge, bool& dir) const {
+ int nid = edge._id;
+ if (dir) {
+ if (nid >= _edge_limit) {
+ nid = (nid - _edge_limit) % (_width - 1) +
+ (nid - _edge_limit) / (_width - 1) * _width;
+ if (nid < _node_num - _width) {
+ edge._id = nid;
+ return;
+ }
+ }
+ if (nid % _width > 0) {
+ edge._id = _edge_limit + nid % _width +
+ (nid / _width) * (_width - 1) - 1;
+ dir = false;
+ return;
+ }
+ if (nid >= _width) {
+ edge._id = nid - _width;
+ dir = false;
+ return;
+ }
+ } else {
+ if (nid >= _edge_limit) {
+ nid = (nid - _edge_limit) % (_width - 1) +
+ (nid - _edge_limit) / (_width - 1) * _width + 1;
+ if (nid >= _width) {
+ edge._id = nid - _width;
+ return;
+ }
+ }
+ }
+ edge._id = -1;
+ dir = true;
+ }
+
+ Arc right(Node n) const {
+ if (n._id % _width < _width - 1) {
+ return Arc(((_edge_limit + n._id % _width +
+ (n._id / _width) * (_width - 1)) << 1) | 1);
+ } else {
+ return INVALID;
+ }
+ }
+
+ Arc left(Node n) const {
+ if (n._id % _width > 0) {
+ return Arc((_edge_limit + n._id % _width +
+ (n._id / _width) * (_width - 1) - 1) << 1);
+ } else {
+ return INVALID;
+ }
+ }
+
+ Arc up(Node n) const {
+ if (n._id < _edge_limit) {
+ return Arc((n._id << 1) | 1);
+ } else {
+ return INVALID;
+ }
+ }
+
+ Arc down(Node n) const {
+ if (n._id >= _width) {
+ return Arc((n._id - _width) << 1);
+ } else {
+ return INVALID;
+ }
+ }
+
+ private:
+ int _width, _height;
+ int _node_num, _edge_num;
+ int _edge_limit;
+ };
+
+
+ typedef GraphExtender ExtendedGridGraphBase;
+
+ /// \ingroup graphs
+ ///
+ /// \brief Grid graph class
+ ///
+ /// This class implements a special graph type. The nodes of the
+ /// graph can be indexed by two integer \c (i,j) value where \c i is
+ /// in the \c [0..width()-1] range and j is in the \c
+ /// [0..height()-1] range. Two nodes are connected in the graph if
+ /// the indexes differ exactly on one position and exactly one is
+ /// the difference. The nodes of the graph can be indexed by position
+ /// with the \c operator()() function. The positions of the nodes can be
+ /// get with \c pos(), \c col() and \c row() members. The outgoing
+ /// arcs can be retrieved with the \c right(), \c up(), \c left()
+ /// and \c down() functions, where the bottom-left corner is the
+ /// origin.
+ ///
+ /// \image html grid_graph.png
+ /// \image latex grid_graph.eps "Grid graph" width=\textwidth
+ ///
+ /// A short example about the basic usage:
+ ///\code
+ /// GridGraph graph(rows, cols);
+ /// GridGraph::NodeMap val(graph);
+ /// for (int i = 0; i < graph.width(); ++i) {
+ /// for (int j = 0; j < graph.height(); ++j) {
+ /// val[graph(i, j)] = i + j;
+ /// }
+ /// }
+ ///\endcode
+ ///
+ /// This graph type is fully conform to the \ref concepts::Graph
+ /// "Graph" concept, and it also has an important extra feature
+ /// that its maps are real \ref concepts::ReferenceMap
+ /// "reference map"s.
+ class GridGraph : public ExtendedGridGraphBase {
+ public:
+
+ typedef ExtendedGridGraphBase Parent;
+
+ /// \brief Map to get the indices of the nodes as dim2::Point.
+ ///
+ /// Map to get the indices of the nodes as dim2::Point.
+ class IndexMap {
+ public:
+ /// \brief The key type of the map
+ typedef GridGraph::Node Key;
+ /// \brief The value type of the map
+ typedef dim2::Point Value;
+
+ /// \brief Constructor
+ ///
+ /// Constructor
+ IndexMap(const GridGraph& graph) : _graph(graph) {}
+
+ /// \brief The subscript operator
+ ///
+ /// The subscript operator.
+ Value operator[](Key key) const {
+ return _graph.pos(key);
+ }
+
+ private:
+ const GridGraph& _graph;
+ };
+
+ /// \brief Map to get the column of the nodes.
+ ///
+ /// Map to get the column of the nodes.
+ class ColMap {
+ public:
+ /// \brief The key type of the map
+ typedef GridGraph::Node Key;
+ /// \brief The value type of the map
+ typedef int Value;
+
+ /// \brief Constructor
+ ///
+ /// Constructor
+ ColMap(const GridGraph& graph) : _graph(graph) {}
+
+ /// \brief The subscript operator
+ ///
+ /// The subscript operator.
+ Value operator[](Key key) const {
+ return _graph.col(key);
+ }
+
+ private:
+ const GridGraph& _graph;
+ };
+
+ /// \brief Map to get the row of the nodes.
+ ///
+ /// Map to get the row of the nodes.
+ class RowMap {
+ public:
+ /// \brief The key type of the map
+ typedef GridGraph::Node Key;
+ /// \brief The value type of the map
+ typedef int Value;
+
+ /// \brief Constructor
+ ///
+ /// Constructor
+ RowMap(const GridGraph& graph) : _graph(graph) {}
+
+ /// \brief The subscript operator
+ ///
+ /// The subscript operator.
+ Value operator[](Key key) const {
+ return _graph.row(key);
+ }
+
+ private:
+ const GridGraph& _graph;
+ };
+
+ /// \brief Constructor
+ ///
+ /// Construct a grid graph with given size.
+ GridGraph(int width, int height) { construct(width, height); }
+
+ /// \brief Resize the graph
+ ///
+ /// Resize the graph. The function will fully destroy and rebuild
+ /// the graph. This cause that the maps of the graph will
+ /// reallocated automatically and the previous values will be
+ /// lost.
+ void resize(int width, int height) {
+ Parent::notifier(Arc()).clear();
+ Parent::notifier(Edge()).clear();
+ Parent::notifier(Node()).clear();
+ construct(width, height);
+ Parent::notifier(Node()).build();
+ Parent::notifier(Edge()).build();
+ Parent::notifier(Arc()).build();
+ }
+
+ /// \brief The node on the given position.
+ ///
+ /// Gives back the node on the given position.
+ Node operator()(int i, int j) const {
+ return Parent::operator()(i, j);
+ }
+
+ /// \brief Gives back the column index of the node.
+ ///
+ /// Gives back the column index of the node.
+ int col(Node n) const {
+ return Parent::col(n);
+ }
+
+ /// \brief Gives back the row index of the node.
+ ///
+ /// Gives back the row index of the node.
+ int row(Node n) const {
+ return Parent::row(n);
+ }
+
+ /// \brief Gives back the position of the node.
+ ///
+ /// Gives back the position of the node, ie. the (col,row) pair.
+ dim2::Point pos(Node n) const {
+ return Parent::pos(n);
+ }
+
+ /// \brief Gives back the number of the columns.
+ ///
+ /// Gives back the number of the columns.
+ int width() const {
+ return Parent::width();
+ }
+
+ /// \brief Gives back the number of the rows.
+ ///
+ /// Gives back the number of the rows.
+ int height() const {
+ return Parent::height();
+ }
+
+ /// \brief Gives back the arc goes right from the node.
+ ///
+ /// Gives back the arc goes right from the node. If there is not
+ /// outgoing arc then it gives back INVALID.
+ Arc right(Node n) const {
+ return Parent::right(n);
+ }
+
+ /// \brief Gives back the arc goes left from the node.
+ ///
+ /// Gives back the arc goes left from the node. If there is not
+ /// outgoing arc then it gives back INVALID.
+ Arc left(Node n) const {
+ return Parent::left(n);
+ }
+
+ /// \brief Gives back the arc goes up from the node.
+ ///
+ /// Gives back the arc goes up from the node. If there is not
+ /// outgoing arc then it gives back INVALID.
+ Arc up(Node n) const {
+ return Parent::up(n);
+ }
+
+ /// \brief Gives back the arc goes down from the node.
+ ///
+ /// Gives back the arc goes down from the node. If there is not
+ /// outgoing arc then it gives back INVALID.
+ Arc down(Node n) const {
+ return Parent::down(n);
+ }
+
+ /// \brief Index map of the grid graph
+ ///
+ /// Just returns an IndexMap for the grid graph.
+ IndexMap indexMap() const {
+ return IndexMap(*this);
+ }
+
+ /// \brief Row map of the grid graph
+ ///
+ /// Just returns a RowMap for the grid graph.
+ RowMap rowMap() const {
+ return RowMap(*this);
+ }
+
+ /// \brief Column map of the grid graph
+ ///
+ /// Just returns a ColMap for the grid graph.
+ ColMap colMap() const {
+ return ColMap(*this);
+ }
+
+ };
+
+}
+#endif
Index: lemon/hypercube_graph.h
===================================================================
--- lemon/hypercube_graph.h (revision 385)
+++ lemon/hypercube_graph.h (revision 385)
@@ -0,0 +1,440 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef HYPERCUBE_GRAPH_H
+#define HYPERCUBE_GRAPH_H
+
+#include
+#include
+#include
+#include
+
+///\ingroup graphs
+///\file
+///\brief HypercubeGraph class.
+
+namespace lemon {
+
+ class HypercubeGraphBase {
+
+ public:
+
+ typedef HypercubeGraphBase Graph;
+
+ class Node;
+ class Edge;
+ class Arc;
+
+ public:
+
+ HypercubeGraphBase() {}
+
+ protected:
+
+ void construct(int dim) {
+ LEMON_ASSERT(dim >= 1, "The number of dimensions must be at least 1.");
+ _dim = dim;
+ _node_num = 1 << dim;
+ _edge_num = dim * (1 << (dim-1));
+ }
+
+ public:
+
+ typedef True NodeNumTag;
+ typedef True EdgeNumTag;
+ typedef True ArcNumTag;
+
+ int nodeNum() const { return _node_num; }
+ int edgeNum() const { return _edge_num; }
+ int arcNum() const { return 2 * _edge_num; }
+
+ int maxNodeId() const { return _node_num - 1; }
+ int maxEdgeId() const { return _edge_num - 1; }
+ int maxArcId() const { return 2 * _edge_num - 1; }
+
+ static Node nodeFromId(int id) { return Node(id); }
+ static Edge edgeFromId(int id) { return Edge(id); }
+ static Arc arcFromId(int id) { return Arc(id); }
+
+ static int id(Node node) { return node._id; }
+ static int id(Edge edge) { return edge._id; }
+ static int id(Arc arc) { return arc._id; }
+
+ Node u(Edge edge) const {
+ int base = edge._id & ((1 << (_dim-1)) - 1);
+ int k = edge._id >> (_dim-1);
+ return ((base >> k) << (k+1)) | (base & ((1 << k) - 1));
+ }
+
+ Node v(Edge edge) const {
+ int base = edge._id & ((1 << (_dim-1)) - 1);
+ int k = edge._id >> (_dim-1);
+ return ((base >> k) << (k+1)) | (base & ((1 << k) - 1)) | (1 << k);
+ }
+
+ Node source(Arc arc) const {
+ return (arc._id & 1) == 1 ? u(arc) : v(arc);
+ }
+
+ Node target(Arc arc) const {
+ return (arc._id & 1) == 1 ? v(arc) : u(arc);
+ }
+
+ typedef True FindEdgeTag;
+ typedef True FindArcTag;
+
+ Edge findEdge(Node u, Node v, Edge prev = INVALID) const {
+ if (prev != INVALID) return INVALID;
+ int d = u._id ^ v._id;
+ int k = 0;
+ if (d == 0) return INVALID;
+ for ( ; (d & 1) == 0; d >>= 1) ++k;
+ if (d >> 1 != 0) return INVALID;
+ return (k << (_dim-1)) | ((u._id >> (k+1)) << k) |
+ (u._id & ((1 << k) - 1));
+ }
+
+ Arc findArc(Node u, Node v, Arc prev = INVALID) const {
+ Edge edge = findEdge(u, v, prev);
+ if (edge == INVALID) return INVALID;
+ int k = edge._id >> (_dim-1);
+ return ((u._id >> k) & 1) == 1 ? edge._id << 1 : (edge._id << 1) | 1;
+ }
+
+ class Node {
+ friend class HypercubeGraphBase;
+
+ protected:
+ int _id;
+ Node(int id) : _id(id) {}
+ public:
+ Node() {}
+ Node (Invalid) : _id(-1) {}
+ bool operator==(const Node node) const {return _id == node._id;}
+ bool operator!=(const Node node) const {return _id != node._id;}
+ bool operator<(const Node node) const {return _id < node._id;}
+ };
+
+ class Edge {
+ friend class HypercubeGraphBase;
+ friend class Arc;
+
+ protected:
+ int _id;
+
+ Edge(int id) : _id(id) {}
+
+ public:
+ Edge() {}
+ Edge (Invalid) : _id(-1) {}
+ bool operator==(const Edge edge) const {return _id == edge._id;}
+ bool operator!=(const Edge edge) const {return _id != edge._id;}
+ bool operator<(const Edge edge) const {return _id < edge._id;}
+ };
+
+ class Arc {
+ friend class HypercubeGraphBase;
+
+ protected:
+ int _id;
+
+ Arc(int id) : _id(id) {}
+
+ public:
+ Arc() {}
+ Arc (Invalid) : _id(-1) {}
+ operator Edge() const { return _id != -1 ? Edge(_id >> 1) : INVALID; }
+ bool operator==(const Arc arc) const {return _id == arc._id;}
+ bool operator!=(const Arc arc) const {return _id != arc._id;}
+ bool operator<(const Arc arc) const {return _id < arc._id;}
+ };
+
+ void first(Node& node) const {
+ node._id = _node_num - 1;
+ }
+
+ static void next(Node& node) {
+ --node._id;
+ }
+
+ void first(Edge& edge) const {
+ edge._id = _edge_num - 1;
+ }
+
+ static void next(Edge& edge) {
+ --edge._id;
+ }
+
+ void first(Arc& arc) const {
+ arc._id = 2 * _edge_num - 1;
+ }
+
+ static void next(Arc& arc) {
+ --arc._id;
+ }
+
+ void firstInc(Edge& edge, bool& dir, const Node& node) const {
+ edge._id = node._id >> 1;
+ dir = (node._id & 1) == 0;
+ }
+
+ void nextInc(Edge& edge, bool& dir) const {
+ Node n = dir ? u(edge) : v(edge);
+ int k = (edge._id >> (_dim-1)) + 1;
+ if (k < _dim) {
+ edge._id = (k << (_dim-1)) |
+ ((n._id >> (k+1)) << k) | (n._id & ((1 << k) - 1));
+ dir = ((n._id >> k) & 1) == 0;
+ } else {
+ edge._id = -1;
+ dir = true;
+ }
+ }
+
+ void firstOut(Arc& arc, const Node& node) const {
+ arc._id = ((node._id >> 1) << 1) | (~node._id & 1);
+ }
+
+ void nextOut(Arc& arc) const {
+ Node n = (arc._id & 1) == 1 ? u(arc) : v(arc);
+ int k = (arc._id >> _dim) + 1;
+ if (k < _dim) {
+ arc._id = (k << (_dim-1)) |
+ ((n._id >> (k+1)) << k) | (n._id & ((1 << k) - 1));
+ arc._id = (arc._id << 1) | (~(n._id >> k) & 1);
+ } else {
+ arc._id = -1;
+ }
+ }
+
+ void firstIn(Arc& arc, const Node& node) const {
+ arc._id = ((node._id >> 1) << 1) | (node._id & 1);
+ }
+
+ void nextIn(Arc& arc) const {
+ Node n = (arc._id & 1) == 1 ? v(arc) : u(arc);
+ int k = (arc._id >> _dim) + 1;
+ if (k < _dim) {
+ arc._id = (k << (_dim-1)) |
+ ((n._id >> (k+1)) << k) | (n._id & ((1 << k) - 1));
+ arc._id = (arc._id << 1) | ((n._id >> k) & 1);
+ } else {
+ arc._id = -1;
+ }
+ }
+
+ static bool direction(Arc arc) {
+ return (arc._id & 1) == 1;
+ }
+
+ static Arc direct(Edge edge, bool dir) {
+ return Arc((edge._id << 1) | (dir ? 1 : 0));
+ }
+
+ int dimension() const {
+ return _dim;
+ }
+
+ bool projection(Node node, int n) const {
+ return static_cast(node._id & (1 << n));
+ }
+
+ int dimension(Edge edge) const {
+ return edge._id >> (_dim-1);
+ }
+
+ int dimension(Arc arc) const {
+ return arc._id >> _dim;
+ }
+
+ int index(Node node) const {
+ return node._id;
+ }
+
+ Node operator()(int ix) const {
+ return Node(ix);
+ }
+
+ private:
+ int _dim;
+ int _node_num, _edge_num;
+ };
+
+
+ typedef GraphExtender ExtendedHypercubeGraphBase;
+
+ /// \ingroup graphs
+ ///
+ /// \brief Hypercube graph class
+ ///
+ /// This class implements a special graph type. The nodes of the graph
+ /// are indiced with integers with at most \c dim binary digits.
+ /// Two nodes are connected in the graph if and only if their indices
+ /// differ only on one position in the binary form.
+ ///
+ /// \note The type of the indices is chosen to \c int for efficiency
+ /// reasons. Thus the maximum dimension of this implementation is 26
+ /// (assuming that the size of \c int is 32 bit).
+ ///
+ /// This graph type is fully conform to the \ref concepts::Graph
+ /// "Graph" concept, and it also has an important extra feature
+ /// that its maps are real \ref concepts::ReferenceMap
+ /// "reference map"s.
+ class HypercubeGraph : public ExtendedHypercubeGraphBase {
+ public:
+
+ typedef ExtendedHypercubeGraphBase Parent;
+
+ /// \brief Constructs a hypercube graph with \c dim dimensions.
+ ///
+ /// Constructs a hypercube graph with \c dim dimensions.
+ HypercubeGraph(int dim) { construct(dim); }
+
+ /// \brief The number of dimensions.
+ ///
+ /// Gives back the number of dimensions.
+ int dimension() const {
+ return Parent::dimension();
+ }
+
+ /// \brief Returns \c true if the n'th bit of the node is one.
+ ///
+ /// Returns \c true if the n'th bit of the node is one.
+ bool projection(Node node, int n) const {
+ return Parent::projection(node, n);
+ }
+
+ /// \brief The dimension id of an edge.
+ ///
+ /// Gives back the dimension id of the given edge.
+ /// It is in the [0..dim-1] range.
+ int dimension(Edge edge) const {
+ return Parent::dimension(edge);
+ }
+
+ /// \brief The dimension id of an arc.
+ ///
+ /// Gives back the dimension id of the given arc.
+ /// It is in the [0..dim-1] range.
+ int dimension(Arc arc) const {
+ return Parent::dimension(arc);
+ }
+
+ /// \brief The index of a node.
+ ///
+ /// Gives back the index of the given node.
+ /// The lower bits of the integer describes the node.
+ int index(Node node) const {
+ return Parent::index(node);
+ }
+
+ /// \brief Gives back a node by its index.
+ ///
+ /// Gives back a node by its index.
+ Node operator()(int ix) const {
+ return Parent::operator()(ix);
+ }
+
+ /// \brief Number of nodes.
+ int nodeNum() const { return Parent::nodeNum(); }
+ /// \brief Number of edges.
+ int edgeNum() const { return Parent::edgeNum(); }
+ /// \brief Number of arcs.
+ int arcNum() const { return Parent::arcNum(); }
+
+ /// \brief Linear combination map.
+ ///
+ /// This map makes possible to give back a linear combination
+ /// for each node. It works like the \c std::accumulate function,
+ /// so it accumulates the \c bf binary function with the \c fv first
+ /// value. The map accumulates only on that positions (dimensions)
+ /// where the index of the node is one. The values that have to be
+ /// accumulated should be given by the \c begin and \c end iterators
+ /// and the length of this range should be equal to the dimension
+ /// number of the graph.
+ ///
+ ///\code
+ /// const int DIM = 3;
+ /// HypercubeGraph graph(DIM);
+ /// dim2::Point base[DIM];
+ /// for (int k = 0; k < DIM; ++k) {
+ /// base[k].x = rnd();
+ /// base[k].y = rnd();
+ /// }
+ /// HypercubeGraph::HyperMap >
+ /// pos(graph, base, base + DIM, dim2::Point(0.0, 0.0));
+ ///\endcode
+ ///
+ /// \see HypercubeGraph
+ template >
+ class HyperMap {
+ public:
+
+ /// \brief The key type of the map
+ typedef Node Key;
+ /// \brief The value type of the map
+ typedef T Value;
+
+ /// \brief Constructor for HyperMap.
+ ///
+ /// Construct a HyperMap for the given graph. The values that have
+ /// to be accumulated should be given by the \c begin and \c end
+ /// iterators and the length of this range should be equal to the
+ /// dimension number of the graph.
+ ///
+ /// This map accumulates the \c bf binary function with the \c fv
+ /// first value on that positions (dimensions) where the index of
+ /// the node is one.
+ template
+ HyperMap(const Graph& graph, It begin, It end,
+ T fv = 0, const BF& bf = BF())
+ : _graph(graph), _values(begin, end), _first_value(fv), _bin_func(bf)
+ {
+ LEMON_ASSERT(_values.size() == graph.dimension(),
+ "Wrong size of range");
+ }
+
+ /// \brief The partial accumulated value.
+ ///
+ /// Gives back the partial accumulated value.
+ Value operator[](const Key& k) const {
+ Value val = _first_value;
+ int id = _graph.index(k);
+ int n = 0;
+ while (id != 0) {
+ if (id & 1) {
+ val = _bin_func(val, _values[n]);
+ }
+ id >>= 1;
+ ++n;
+ }
+ return val;
+ }
+
+ private:
+ const Graph& _graph;
+ std::vector _values;
+ T _first_value;
+ BF _bin_func;
+ };
+
+ };
+
+}
+
+#endif
Index: lemon/list_graph.h
===================================================================
--- lemon/list_graph.h (revision 313)
+++ lemon/list_graph.h (revision 341)
@@ -841,6 +841,6 @@
public:
- operator Edge() const {
- return id != -1 ? edgeFromId(id / 2) : INVALID;
+ operator Edge() const {
+ return id != -1 ? edgeFromId(id / 2) : INVALID;
}
Index: lemon/max_matching.h
===================================================================
--- lemon/max_matching.h (revision 342)
+++ lemon/max_matching.h (revision 342)
@@ -0,0 +1,3104 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_MAX_MATCHING_H
+#define LEMON_MAX_MATCHING_H
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+///\ingroup matching
+///\file
+///\brief Maximum matching algorithms in general graphs.
+
+namespace lemon {
+
+ /// \ingroup matching
+ ///
+ /// \brief Edmonds' alternating forest maximum matching algorithm.
+ ///
+ /// This class implements Edmonds' alternating forest matching
+ /// algorithm. The algorithm can be started from an arbitrary initial
+ /// matching (the default is the empty one)
+ ///
+ /// The dual solution of the problem is a map of the nodes to
+ /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
+ /// MATCHED/C showing the Gallai-Edmonds decomposition of the
+ /// graph. The nodes in \c EVEN/D induce a graph with
+ /// factor-critical components, the nodes in \c ODD/A form the
+ /// barrier, and the nodes in \c MATCHED/C induce a graph having a
+ /// perfect matching. The number of the factor-critical components
+ /// minus the number of barrier nodes is a lower bound on the
+ /// unmatched nodes, and the matching is optimal if and only if this bound is
+ /// tight. This decomposition can be attained by calling \c
+ /// decomposition() after running the algorithm.
+ ///
+ /// \param _Graph The graph type the algorithm runs on.
+ template
+ class MaxMatching {
+ public:
+
+ typedef _Graph Graph;
+ typedef typename Graph::template NodeMap
+ MatchingMap;
+
+ ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
+ ///
+ ///Indicates the Gallai-Edmonds decomposition of the graph. The
+ ///nodes with Status \c EVEN/D induce a graph with factor-critical
+ ///components, the nodes in \c ODD/A form the canonical barrier,
+ ///and the nodes in \c MATCHED/C induce a graph having a perfect
+ ///matching.
+ enum Status {
+ EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
+ };
+
+ typedef typename Graph::template NodeMap StatusMap;
+
+ private:
+
+ TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+ typedef UnionFindEnum BlossomSet;
+ typedef ExtendFindEnum TreeSet;
+ typedef RangeMap NodeIntMap;
+ typedef MatchingMap EarMap;
+ typedef std::vector NodeQueue;
+
+ const Graph& _graph;
+ MatchingMap* _matching;
+ StatusMap* _status;
+
+ EarMap* _ear;
+
+ IntNodeMap* _blossom_set_index;
+ BlossomSet* _blossom_set;
+ NodeIntMap* _blossom_rep;
+
+ IntNodeMap* _tree_set_index;
+ TreeSet* _tree_set;
+
+ NodeQueue _node_queue;
+ int _process, _postpone, _last;
+
+ int _node_num;
+
+ private:
+
+ void createStructures() {
+ _node_num = countNodes(_graph);
+ if (!_matching) {
+ _matching = new MatchingMap(_graph);
+ }
+ if (!_status) {
+ _status = new StatusMap(_graph);
+ }
+ if (!_ear) {
+ _ear = new EarMap(_graph);
+ }
+ if (!_blossom_set) {
+ _blossom_set_index = new IntNodeMap(_graph);
+ _blossom_set = new BlossomSet(*_blossom_set_index);
+ }
+ if (!_blossom_rep) {
+ _blossom_rep = new NodeIntMap(_node_num);
+ }
+ if (!_tree_set) {
+ _tree_set_index = new IntNodeMap(_graph);
+ _tree_set = new TreeSet(*_tree_set_index);
+ }
+ _node_queue.resize(_node_num);
+ }
+
+ void destroyStructures() {
+ if (_matching) {
+ delete _matching;
+ }
+ if (_status) {
+ delete _status;
+ }
+ if (_ear) {
+ delete _ear;
+ }
+ if (_blossom_set) {
+ delete _blossom_set;
+ delete _blossom_set_index;
+ }
+ if (_blossom_rep) {
+ delete _blossom_rep;
+ }
+ if (_tree_set) {
+ delete _tree_set_index;
+ delete _tree_set;
+ }
+ }
+
+ void processDense(const Node& n) {
+ _process = _postpone = _last = 0;
+ _node_queue[_last++] = n;
+
+ while (_process != _last) {
+ Node u = _node_queue[_process++];
+ for (OutArcIt a(_graph, u); a != INVALID; ++a) {
+ Node v = _graph.target(a);
+ if ((*_status)[v] == MATCHED) {
+ extendOnArc(a);
+ } else if ((*_status)[v] == UNMATCHED) {
+ augmentOnArc(a);
+ return;
+ }
+ }
+ }
+
+ while (_postpone != _last) {
+ Node u = _node_queue[_postpone++];
+
+ for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
+ Node v = _graph.target(a);
+
+ if ((*_status)[v] == EVEN) {
+ if (_blossom_set->find(u) != _blossom_set->find(v)) {
+ shrinkOnEdge(a);
+ }
+ }
+
+ while (_process != _last) {
+ Node w = _node_queue[_process++];
+ for (OutArcIt b(_graph, w); b != INVALID; ++b) {
+ Node x = _graph.target(b);
+ if ((*_status)[x] == MATCHED) {
+ extendOnArc(b);
+ } else if ((*_status)[x] == UNMATCHED) {
+ augmentOnArc(b);
+ return;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void processSparse(const Node& n) {
+ _process = _last = 0;
+ _node_queue[_last++] = n;
+ while (_process != _last) {
+ Node u = _node_queue[_process++];
+ for (OutArcIt a(_graph, u); a != INVALID; ++a) {
+ Node v = _graph.target(a);
+
+ if ((*_status)[v] == EVEN) {
+ if (_blossom_set->find(u) != _blossom_set->find(v)) {
+ shrinkOnEdge(a);
+ }
+ } else if ((*_status)[v] == MATCHED) {
+ extendOnArc(a);
+ } else if ((*_status)[v] == UNMATCHED) {
+ augmentOnArc(a);
+ return;
+ }
+ }
+ }
+ }
+
+ void shrinkOnEdge(const Edge& e) {
+ Node nca = INVALID;
+
+ {
+ std::set left_set, right_set;
+
+ Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
+ left_set.insert(left);
+
+ Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
+ right_set.insert(right);
+
+ while (true) {
+ if ((*_matching)[left] == INVALID) break;
+ left = _graph.target((*_matching)[left]);
+ left = (*_blossom_rep)[_blossom_set->
+ find(_graph.target((*_ear)[left]))];
+ if (right_set.find(left) != right_set.end()) {
+ nca = left;
+ break;
+ }
+ left_set.insert(left);
+
+ if ((*_matching)[right] == INVALID) break;
+ right = _graph.target((*_matching)[right]);
+ right = (*_blossom_rep)[_blossom_set->
+ find(_graph.target((*_ear)[right]))];
+ if (left_set.find(right) != left_set.end()) {
+ nca = right;
+ break;
+ }
+ right_set.insert(right);
+ }
+
+ if (nca == INVALID) {
+ if ((*_matching)[left] == INVALID) {
+ nca = right;
+ while (left_set.find(nca) == left_set.end()) {
+ nca = _graph.target((*_matching)[nca]);
+ nca =(*_blossom_rep)[_blossom_set->
+ find(_graph.target((*_ear)[nca]))];
+ }
+ } else {
+ nca = left;
+ while (right_set.find(nca) == right_set.end()) {
+ nca = _graph.target((*_matching)[nca]);
+ nca = (*_blossom_rep)[_blossom_set->
+ find(_graph.target((*_ear)[nca]))];
+ }
+ }
+ }
+ }
+
+ {
+
+ Node node = _graph.u(e);
+ Arc arc = _graph.direct(e, true);
+ Node base = (*_blossom_rep)[_blossom_set->find(node)];
+
+ while (base != nca) {
+ _ear->set(node, arc);
+
+ Node n = node;
+ while (n != base) {
+ n = _graph.target((*_matching)[n]);
+ Arc a = (*_ear)[n];
+ n = _graph.target(a);
+ _ear->set(n, _graph.oppositeArc(a));
+ }
+ node = _graph.target((*_matching)[base]);
+ _tree_set->erase(base);
+ _tree_set->erase(node);
+ _blossom_set->insert(node, _blossom_set->find(base));
+ _status->set(node, EVEN);
+ _node_queue[_last++] = node;
+ arc = _graph.oppositeArc((*_ear)[node]);
+ node = _graph.target((*_ear)[node]);
+ base = (*_blossom_rep)[_blossom_set->find(node)];
+ _blossom_set->join(_graph.target(arc), base);
+ }
+ }
+
+ _blossom_rep->set(_blossom_set->find(nca), nca);
+
+ {
+
+ Node node = _graph.v(e);
+ Arc arc = _graph.direct(e, false);
+ Node base = (*_blossom_rep)[_blossom_set->find(node)];
+
+ while (base != nca) {
+ _ear->set(node, arc);
+
+ Node n = node;
+ while (n != base) {
+ n = _graph.target((*_matching)[n]);
+ Arc a = (*_ear)[n];
+ n = _graph.target(a);
+ _ear->set(n, _graph.oppositeArc(a));
+ }
+ node = _graph.target((*_matching)[base]);
+ _tree_set->erase(base);
+ _tree_set->erase(node);
+ _blossom_set->insert(node, _blossom_set->find(base));
+ _status->set(node, EVEN);
+ _node_queue[_last++] = node;
+ arc = _graph.oppositeArc((*_ear)[node]);
+ node = _graph.target((*_ear)[node]);
+ base = (*_blossom_rep)[_blossom_set->find(node)];
+ _blossom_set->join(_graph.target(arc), base);
+ }
+ }
+
+ _blossom_rep->set(_blossom_set->find(nca), nca);
+ }
+
+
+
+ void extendOnArc(const Arc& a) {
+ Node base = _graph.source(a);
+ Node odd = _graph.target(a);
+
+ _ear->set(odd, _graph.oppositeArc(a));
+ Node even = _graph.target((*_matching)[odd]);
+ _blossom_rep->set(_blossom_set->insert(even), even);
+ _status->set(odd, ODD);
+ _status->set(even, EVEN);
+ int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
+ _tree_set->insert(odd, tree);
+ _tree_set->insert(even, tree);
+ _node_queue[_last++] = even;
+
+ }
+
+ void augmentOnArc(const Arc& a) {
+ Node even = _graph.source(a);
+ Node odd = _graph.target(a);
+
+ int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
+
+ _matching->set(odd, _graph.oppositeArc(a));
+ _status->set(odd, MATCHED);
+
+ Arc arc = (*_matching)[even];
+ _matching->set(even, a);
+
+ while (arc != INVALID) {
+ odd = _graph.target(arc);
+ arc = (*_ear)[odd];
+ even = _graph.target(arc);
+ _matching->set(odd, arc);
+ arc = (*_matching)[even];
+ _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
+ }
+
+ for (typename TreeSet::ItemIt it(*_tree_set, tree);
+ it != INVALID; ++it) {
+ if ((*_status)[it] == ODD) {
+ _status->set(it, MATCHED);
+ } else {
+ int blossom = _blossom_set->find(it);
+ for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
+ jt != INVALID; ++jt) {
+ _status->set(jt, MATCHED);
+ }
+ _blossom_set->eraseClass(blossom);
+ }
+ }
+ _tree_set->eraseClass(tree);
+
+ }
+
+ public:
+
+ /// \brief Constructor
+ ///
+ /// Constructor.
+ MaxMatching(const Graph& graph)
+ : _graph(graph), _matching(0), _status(0), _ear(0),
+ _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
+ _tree_set_index(0), _tree_set(0) {}
+
+ ~MaxMatching() {
+ destroyStructures();
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to use the
+ /// \c run() member function.
+ /// \n
+
+ /// If you need better control on the execution, you must call
+ /// \ref init(), \ref greedyInit() or \ref matchingInit()
+ /// functions first, then you can start the algorithm with the \ref
+ /// startParse() or startDense() functions.
+
+ ///@{
+
+ /// \brief Sets the actual matching to the empty matching.
+ ///
+ /// Sets the actual matching to the empty matching.
+ ///
+ void init() {
+ createStructures();
+ for(NodeIt n(_graph); n != INVALID; ++n) {
+ _matching->set(n, INVALID);
+ _status->set(n, UNMATCHED);
+ }
+ }
+
+ ///\brief Finds an initial matching in a greedy way
+ ///
+ ///It finds an initial matching in a greedy way.
+ void greedyInit() {
+ createStructures();
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _matching->set(n, INVALID);
+ _status->set(n, UNMATCHED);
+ }
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_matching)[n] == INVALID) {
+ for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
+ Node v = _graph.target(a);
+ if ((*_matching)[v] == INVALID && v != n) {
+ _matching->set(n, a);
+ _status->set(n, MATCHED);
+ _matching->set(v, _graph.oppositeArc(a));
+ _status->set(v, MATCHED);
+ break;
+ }
+ }
+ }
+ }
+ }
+
+
+ /// \brief Initialize the matching from a map containing.
+ ///
+ /// Initialize the matching from a \c bool valued \c Edge map. This
+ /// map must have the property that there are no two incident edges
+ /// with true value, ie. it contains a matching.
+ /// \return %True if the map contains a matching.
+ template
+ bool matchingInit(const MatchingMap& matching) {
+ createStructures();
+
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _matching->set(n, INVALID);
+ _status->set(n, UNMATCHED);
+ }
+ for(EdgeIt e(_graph); e!=INVALID; ++e) {
+ if (matching[e]) {
+
+ Node u = _graph.u(e);
+ if ((*_matching)[u] != INVALID) return false;
+ _matching->set(u, _graph.direct(e, true));
+ _status->set(u, MATCHED);
+
+ Node v = _graph.v(e);
+ if ((*_matching)[v] != INVALID) return false;
+ _matching->set(v, _graph.direct(e, false));
+ _status->set(v, MATCHED);
+ }
+ }
+ return true;
+ }
+
+ /// \brief Starts Edmonds' algorithm
+ ///
+ /// If runs the original Edmonds' algorithm.
+ void startSparse() {
+ for(NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_status)[n] == UNMATCHED) {
+ (*_blossom_rep)[_blossom_set->insert(n)] = n;
+ _tree_set->insert(n);
+ _status->set(n, EVEN);
+ processSparse(n);
+ }
+ }
+ }
+
+ /// \brief Starts Edmonds' algorithm.
+ ///
+ /// It runs Edmonds' algorithm with a heuristic of postponing
+ /// shrinks, therefore resulting in a faster algorithm for dense graphs.
+ void startDense() {
+ for(NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_status)[n] == UNMATCHED) {
+ (*_blossom_rep)[_blossom_set->insert(n)] = n;
+ _tree_set->insert(n);
+ _status->set(n, EVEN);
+ processDense(n);
+ }
+ }
+ }
+
+
+ /// \brief Runs Edmonds' algorithm
+ ///
+ /// Runs Edmonds' algorithm for sparse graphs (m<2*n)
+ /// or Edmonds' algorithm with a heuristic of
+ /// postponing shrinks for dense graphs.
+ void run() {
+ if (countEdges(_graph) < 2 * countNodes(_graph)) {
+ greedyInit();
+ startSparse();
+ } else {
+ init();
+ startDense();
+ }
+ }
+
+ /// @}
+
+ /// \name Primal solution
+ /// Functions to get the primal solution, ie. the matching.
+
+ /// @{
+
+ ///\brief Returns the size of the current matching.
+ ///
+ ///Returns the size of the current matching. After \ref
+ ///run() it returns the size of the maximum matching in the graph.
+ int matchingSize() const {
+ int size = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_matching)[n] != INVALID) {
+ ++size;
+ }
+ }
+ return size / 2;
+ }
+
+ /// \brief Returns true when the edge is in the matching.
+ ///
+ /// Returns true when the edge is in the matching.
+ bool matching(const Edge& edge) const {
+ return edge == (*_matching)[_graph.u(edge)];
+ }
+
+ /// \brief Returns the matching edge incident to the given node.
+ ///
+ /// Returns the matching edge of a \c node in the actual matching or
+ /// INVALID if the \c node is not covered by the actual matching.
+ Arc matching(const Node& n) const {
+ return (*_matching)[n];
+ }
+
+ ///\brief Returns the mate of a node in the actual matching.
+ ///
+ ///Returns the mate of a \c node in the actual matching or
+ ///INVALID if the \c node is not covered by the actual matching.
+ Node mate(const Node& n) const {
+ return (*_matching)[n] != INVALID ?
+ _graph.target((*_matching)[n]) : INVALID;
+ }
+
+ /// @}
+
+ /// \name Dual solution
+ /// Functions to get the dual solution, ie. the decomposition.
+
+ /// @{
+
+ /// \brief Returns the class of the node in the Edmonds-Gallai
+ /// decomposition.
+ ///
+ /// Returns the class of the node in the Edmonds-Gallai
+ /// decomposition.
+ Status decomposition(const Node& n) const {
+ return (*_status)[n];
+ }
+
+ /// \brief Returns true when the node is in the barrier.
+ ///
+ /// Returns true when the node is in the barrier.
+ bool barrier(const Node& n) const {
+ return (*_status)[n] == ODD;
+ }
+
+ /// @}
+
+ };
+
+ /// \ingroup matching
+ ///
+ /// \brief Weighted matching in general graphs
+ ///
+ /// This class provides an efficient implementation of Edmond's
+ /// maximum weighted matching algorithm. The implementation is based
+ /// on extensive use of priority queues and provides
+ /// \f$O(nm\log(n))\f$ time complexity.
+ ///
+ /// The maximum weighted matching problem is to find undirected
+ /// edges in the graph with maximum overall weight and no two of
+ /// them shares their ends. The problem can be formulated with the
+ /// following linear program.
+ /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
+ /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
+ \quad \forall B\in\mathcal{O}\f] */
+ /// \f[x_e \ge 0\quad \forall e\in E\f]
+ /// \f[\max \sum_{e\in E}x_ew_e\f]
+ /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+ /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
+ /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
+ /// subsets of the nodes.
+ ///
+ /// The algorithm calculates an optimal matching and a proof of the
+ /// optimality. The solution of the dual problem can be used to check
+ /// the result of the algorithm. The dual linear problem is the
+ /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
+ z_B \ge w_{uv} \quad \forall uv\in E\f] */
+ /// \f[y_u \ge 0 \quad \forall u \in V\f]
+ /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
+ /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
+ \frac{\vert B \vert - 1}{2}z_B\f] */
+ ///
+ /// The algorithm can be executed with \c run() or the \c init() and
+ /// then the \c start() member functions. After it the matching can
+ /// be asked with \c matching() or mate() functions. The dual
+ /// solution can be get with \c nodeValue(), \c blossomNum() and \c
+ /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
+ /// "BlossomIt" nested class, which is able to iterate on the nodes
+ /// of a blossom. If the value type is integral then the dual
+ /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
+ template >
+ class MaxWeightedMatching {
+ public:
+
+ typedef _Graph Graph;
+ typedef _WeightMap WeightMap;
+ typedef typename WeightMap::Value Value;
+
+ /// \brief Scaling factor for dual solution
+ ///
+ /// Scaling factor for dual solution, it is equal to 4 or 1
+ /// according to the value type.
+ static const int dualScale =
+ std::numeric_limits::is_integer ? 4 : 1;
+
+ typedef typename Graph::template NodeMap
+ MatchingMap;
+
+ private:
+
+ TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+ typedef typename Graph::template NodeMap NodePotential;
+ typedef std::vector BlossomNodeList;
+
+ struct BlossomVariable {
+ int begin, end;
+ Value value;
+
+ BlossomVariable(int _begin, int _end, Value _value)
+ : begin(_begin), end(_end), value(_value) {}
+
+ };
+
+ typedef std::vector BlossomPotential;
+
+ const Graph& _graph;
+ const WeightMap& _weight;
+
+ MatchingMap* _matching;
+
+ NodePotential* _node_potential;
+
+ BlossomPotential _blossom_potential;
+ BlossomNodeList _blossom_node_list;
+
+ int _node_num;
+ int _blossom_num;
+
+ typedef RangeMap IntIntMap;
+
+ enum Status {
+ EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
+ };
+
+ typedef HeapUnionFind BlossomSet;
+ struct BlossomData {
+ int tree;
+ Status status;
+ Arc pred, next;
+ Value pot, offset;
+ Node base;
+ };
+
+ IntNodeMap *_blossom_index;
+ BlossomSet *_blossom_set;
+ RangeMap* _blossom_data;
+
+ IntNodeMap *_node_index;
+ IntArcMap *_node_heap_index;
+
+ struct NodeData {
+
+ NodeData(IntArcMap& node_heap_index)
+ : heap(node_heap_index) {}
+
+ int blossom;
+ Value pot;
+ BinHeap heap;
+ std::map heap_index;
+
+ int tree;
+ };
+
+ RangeMap* _node_data;
+
+ typedef ExtendFindEnum TreeSet;
+
+ IntIntMap *_tree_set_index;
+ TreeSet *_tree_set;
+
+ IntNodeMap *_delta1_index;
+ BinHeap *_delta1;
+
+ IntIntMap *_delta2_index;
+ BinHeap *_delta2;
+
+ IntEdgeMap *_delta3_index;
+ BinHeap *_delta3;
+
+ IntIntMap *_delta4_index;
+ BinHeap *_delta4;
+
+ Value _delta_sum;
+
+ void createStructures() {
+ _node_num = countNodes(_graph);
+ _blossom_num = _node_num * 3 / 2;
+
+ if (!_matching) {
+ _matching = new MatchingMap(_graph);
+ }
+ if (!_node_potential) {
+ _node_potential = new NodePotential(_graph);
+ }
+ if (!_blossom_set) {
+ _blossom_index = new IntNodeMap(_graph);
+ _blossom_set = new BlossomSet(*_blossom_index);
+ _blossom_data = new RangeMap(_blossom_num);
+ }
+
+ if (!_node_index) {
+ _node_index = new IntNodeMap(_graph);
+ _node_heap_index = new IntArcMap(_graph);
+ _node_data = new RangeMap(_node_num,
+ NodeData(*_node_heap_index));
+ }
+
+ if (!_tree_set) {
+ _tree_set_index = new IntIntMap(_blossom_num);
+ _tree_set = new TreeSet(*_tree_set_index);
+ }
+ if (!_delta1) {
+ _delta1_index = new IntNodeMap(_graph);
+ _delta1 = new BinHeap(*_delta1_index);
+ }
+ if (!_delta2) {
+ _delta2_index = new IntIntMap(_blossom_num);
+ _delta2 = new BinHeap(*_delta2_index);
+ }
+ if (!_delta3) {
+ _delta3_index = new IntEdgeMap(_graph);
+ _delta3 = new BinHeap(*_delta3_index);
+ }
+ if (!_delta4) {
+ _delta4_index = new IntIntMap(_blossom_num);
+ _delta4 = new BinHeap(*_delta4_index);
+ }
+ }
+
+ void destroyStructures() {
+ _node_num = countNodes(_graph);
+ _blossom_num = _node_num * 3 / 2;
+
+ if (_matching) {
+ delete _matching;
+ }
+ if (_node_potential) {
+ delete _node_potential;
+ }
+ if (_blossom_set) {
+ delete _blossom_index;
+ delete _blossom_set;
+ delete _blossom_data;
+ }
+
+ if (_node_index) {
+ delete _node_index;
+ delete _node_heap_index;
+ delete _node_data;
+ }
+
+ if (_tree_set) {
+ delete _tree_set_index;
+ delete _tree_set;
+ }
+ if (_delta1) {
+ delete _delta1_index;
+ delete _delta1;
+ }
+ if (_delta2) {
+ delete _delta2_index;
+ delete _delta2;
+ }
+ if (_delta3) {
+ delete _delta3_index;
+ delete _delta3;
+ }
+ if (_delta4) {
+ delete _delta4_index;
+ delete _delta4;
+ }
+ }
+
+ void matchedToEven(int blossom, int tree) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+
+ if (!_blossom_set->trivial(blossom)) {
+ (*_blossom_data)[blossom].pot -=
+ 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+
+ _blossom_set->increase(n, std::numeric_limits::max());
+ int ni = (*_node_index)[n];
+
+ (*_node_data)[ni].heap.clear();
+ (*_node_data)[ni].heap_index.clear();
+
+ (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
+
+ _delta1->push(n, (*_node_data)[ni].pot);
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if ((*_blossom_data)[vb].status == EVEN) {
+ if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+ _delta3->push(e, rw / 2);
+ }
+ } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+ if (_delta3->state(e) != _delta3->IN_HEAP) {
+ _delta3->push(e, rw);
+ }
+ } else {
+ typename std::map::iterator it =
+ (*_node_data)[vi].heap_index.find(tree);
+
+ if (it != (*_node_data)[vi].heap_index.end()) {
+ if ((*_node_data)[vi].heap[it->second] > rw) {
+ (*_node_data)[vi].heap.replace(it->second, e);
+ (*_node_data)[vi].heap.decrease(e, rw);
+ it->second = e;
+ }
+ } else {
+ (*_node_data)[vi].heap.push(e, rw);
+ (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+ }
+
+ if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+ _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+ if ((*_blossom_data)[vb].status == MATCHED) {
+ if (_delta2->state(vb) != _delta2->IN_HEAP) {
+ _delta2->push(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset){
+ _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ }
+ }
+ }
+ }
+ }
+ }
+ (*_blossom_data)[blossom].offset = 0;
+ }
+
+ void matchedToOdd(int blossom) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+ (*_blossom_data)[blossom].offset += _delta_sum;
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+
+ void evenToMatched(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ (*_blossom_data)[blossom].pot += 2 * _delta_sum;
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+ (*_node_data)[ni].pot -= _delta_sum;
+
+ _delta1->erase(n);
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if (vb == blossom) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else if ((*_blossom_data)[vb].status == EVEN) {
+
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+
+ int vt = _tree_set->find(vb);
+
+ if (vt != tree) {
+
+ Arc r = _graph.oppositeArc(e);
+
+ typename std::map::iterator it =
+ (*_node_data)[ni].heap_index.find(vt);
+
+ if (it != (*_node_data)[ni].heap_index.end()) {
+ if ((*_node_data)[ni].heap[it->second] > rw) {
+ (*_node_data)[ni].heap.replace(it->second, r);
+ (*_node_data)[ni].heap.decrease(r, rw);
+ it->second = r;
+ }
+ } else {
+ (*_node_data)[ni].heap.push(r, rw);
+ (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+ }
+
+ if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
+ _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+
+ if (_delta2->state(blossom) != _delta2->IN_HEAP) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ } else if ((*_delta2)[blossom] >
+ _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset){
+ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+ }
+
+ } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else {
+
+ typename std::map::iterator it =
+ (*_node_data)[vi].heap_index.find(tree);
+
+ if (it != (*_node_data)[vi].heap_index.end()) {
+ (*_node_data)[vi].heap.erase(it->second);
+ (*_node_data)[vi].heap_index.erase(it);
+ if ((*_node_data)[vi].heap.empty()) {
+ _blossom_set->increase(v, std::numeric_limits::max());
+ } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
+ _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
+ }
+
+ if ((*_blossom_data)[vb].status == MATCHED) {
+ if (_blossom_set->classPrio(vb) ==
+ std::numeric_limits::max()) {
+ _delta2->erase(vb);
+ } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset) {
+ _delta2->increase(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void oddToMatched(int blossom) {
+ (*_blossom_data)[blossom].offset -= _delta_sum;
+
+ if (_blossom_set->classPrio(blossom) !=
+ std::numeric_limits::max()) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ }
+ }
+
+ void oddToEven(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ (*_blossom_data)[blossom].pot -=
+ 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+
+ _blossom_set->increase(n, std::numeric_limits::max());
+
+ (*_node_data)[ni].heap.clear();
+ (*_node_data)[ni].heap_index.clear();
+ (*_node_data)[ni].pot +=
+ 2 * _delta_sum - (*_blossom_data)[blossom].offset;
+
+ _delta1->push(n, (*_node_data)[ni].pot);
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if ((*_blossom_data)[vb].status == EVEN) {
+ if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+ _delta3->push(e, rw / 2);
+ }
+ } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+ if (_delta3->state(e) != _delta3->IN_HEAP) {
+ _delta3->push(e, rw);
+ }
+ } else {
+
+ typename std::map::iterator it =
+ (*_node_data)[vi].heap_index.find(tree);
+
+ if (it != (*_node_data)[vi].heap_index.end()) {
+ if ((*_node_data)[vi].heap[it->second] > rw) {
+ (*_node_data)[vi].heap.replace(it->second, e);
+ (*_node_data)[vi].heap.decrease(e, rw);
+ it->second = e;
+ }
+ } else {
+ (*_node_data)[vi].heap.push(e, rw);
+ (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+ }
+
+ if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+ _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+ if ((*_blossom_data)[vb].status == MATCHED) {
+ if (_delta2->state(vb) != _delta2->IN_HEAP) {
+ _delta2->push(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset) {
+ _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ }
+ }
+ }
+ }
+ }
+ }
+ (*_blossom_data)[blossom].offset = 0;
+ }
+
+
+ void matchedToUnmatched(int blossom) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+
+ _blossom_set->increase(n, std::numeric_limits::max());
+
+ (*_node_data)[ni].heap.clear();
+ (*_node_data)[ni].heap_index.clear();
+
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if ((*_blossom_data)[vb].status == EVEN) {
+ if (_delta3->state(e) != _delta3->IN_HEAP) {
+ _delta3->push(e, rw);
+ }
+ }
+ }
+ }
+ }
+
+ void unmatchedToMatched(int blossom) {
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if (vb == blossom) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else if ((*_blossom_data)[vb].status == EVEN) {
+
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+
+ int vt = _tree_set->find(vb);
+
+ Arc r = _graph.oppositeArc(e);
+
+ typename std::map::iterator it =
+ (*_node_data)[ni].heap_index.find(vt);
+
+ if (it != (*_node_data)[ni].heap_index.end()) {
+ if ((*_node_data)[ni].heap[it->second] > rw) {
+ (*_node_data)[ni].heap.replace(it->second, r);
+ (*_node_data)[ni].heap.decrease(r, rw);
+ it->second = r;
+ }
+ } else {
+ (*_node_data)[ni].heap.push(r, rw);
+ (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+ }
+
+ if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
+ _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+
+ if (_delta2->state(blossom) != _delta2->IN_HEAP) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
+ (*_blossom_data)[blossom].offset){
+ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+
+ } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ }
+ }
+ }
+ }
+
+ void alternatePath(int even, int tree) {
+ int odd;
+
+ evenToMatched(even, tree);
+ (*_blossom_data)[even].status = MATCHED;
+
+ while ((*_blossom_data)[even].pred != INVALID) {
+ odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
+ (*_blossom_data)[odd].status = MATCHED;
+ oddToMatched(odd);
+ (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
+
+ even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
+ (*_blossom_data)[even].status = MATCHED;
+ evenToMatched(even, tree);
+ (*_blossom_data)[even].next =
+ _graph.oppositeArc((*_blossom_data)[odd].pred);
+ }
+
+ }
+
+ void destroyTree(int tree) {
+ for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
+ if ((*_blossom_data)[b].status == EVEN) {
+ (*_blossom_data)[b].status = MATCHED;
+ evenToMatched(b, tree);
+ } else if ((*_blossom_data)[b].status == ODD) {
+ (*_blossom_data)[b].status = MATCHED;
+ oddToMatched(b);
+ }
+ }
+ _tree_set->eraseClass(tree);
+ }
+
+
+ void unmatchNode(const Node& node) {
+ int blossom = _blossom_set->find(node);
+ int tree = _tree_set->find(blossom);
+
+ alternatePath(blossom, tree);
+ destroyTree(tree);
+
+ (*_blossom_data)[blossom].status = UNMATCHED;
+ (*_blossom_data)[blossom].base = node;
+ matchedToUnmatched(blossom);
+ }
+
+
+ void augmentOnEdge(const Edge& edge) {
+
+ int left = _blossom_set->find(_graph.u(edge));
+ int right = _blossom_set->find(_graph.v(edge));
+
+ if ((*_blossom_data)[left].status == EVEN) {
+ int left_tree = _tree_set->find(left);
+ alternatePath(left, left_tree);
+ destroyTree(left_tree);
+ } else {
+ (*_blossom_data)[left].status = MATCHED;
+ unmatchedToMatched(left);
+ }
+
+ if ((*_blossom_data)[right].status == EVEN) {
+ int right_tree = _tree_set->find(right);
+ alternatePath(right, right_tree);
+ destroyTree(right_tree);
+ } else {
+ (*_blossom_data)[right].status = MATCHED;
+ unmatchedToMatched(right);
+ }
+
+ (*_blossom_data)[left].next = _graph.direct(edge, true);
+ (*_blossom_data)[right].next = _graph.direct(edge, false);
+ }
+
+ void extendOnArc(const Arc& arc) {
+ int base = _blossom_set->find(_graph.target(arc));
+ int tree = _tree_set->find(base);
+
+ int odd = _blossom_set->find(_graph.source(arc));
+ _tree_set->insert(odd, tree);
+ (*_blossom_data)[odd].status = ODD;
+ matchedToOdd(odd);
+ (*_blossom_data)[odd].pred = arc;
+
+ int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
+ (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
+ _tree_set->insert(even, tree);
+ (*_blossom_data)[even].status = EVEN;
+ matchedToEven(even, tree);
+ }
+
+ void shrinkOnEdge(const Edge& edge, int tree) {
+ int nca = -1;
+ std::vector left_path, right_path;
+
+ {
+ std::set left_set, right_set;
+ int left = _blossom_set->find(_graph.u(edge));
+ left_path.push_back(left);
+ left_set.insert(left);
+
+ int right = _blossom_set->find(_graph.v(edge));
+ right_path.push_back(right);
+ right_set.insert(right);
+
+ while (true) {
+
+ if ((*_blossom_data)[left].pred == INVALID) break;
+
+ left =
+ _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
+ left_path.push_back(left);
+ left =
+ _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
+ left_path.push_back(left);
+
+ left_set.insert(left);
+
+ if (right_set.find(left) != right_set.end()) {
+ nca = left;
+ break;
+ }
+
+ if ((*_blossom_data)[right].pred == INVALID) break;
+
+ right =
+ _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
+ right_path.push_back(right);
+ right =
+ _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
+ right_path.push_back(right);
+
+ right_set.insert(right);
+
+ if (left_set.find(right) != left_set.end()) {
+ nca = right;
+ break;
+ }
+
+ }
+
+ if (nca == -1) {
+ if ((*_blossom_data)[left].pred == INVALID) {
+ nca = right;
+ while (left_set.find(nca) == left_set.end()) {
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ right_path.push_back(nca);
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ right_path.push_back(nca);
+ }
+ } else {
+ nca = left;
+ while (right_set.find(nca) == right_set.end()) {
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ }
+ }
+ }
+ }
+
+ std::vector subblossoms;
+ Arc prev;
+
+ prev = _graph.direct(edge, true);
+ for (int i = 0; left_path[i] != nca; i += 2) {
+ subblossoms.push_back(left_path[i]);
+ (*_blossom_data)[left_path[i]].next = prev;
+ _tree_set->erase(left_path[i]);
+
+ subblossoms.push_back(left_path[i + 1]);
+ (*_blossom_data)[left_path[i + 1]].status = EVEN;
+ oddToEven(left_path[i + 1], tree);
+ _tree_set->erase(left_path[i + 1]);
+ prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
+ }
+
+ int k = 0;
+ while (right_path[k] != nca) ++k;
+
+ subblossoms.push_back(nca);
+ (*_blossom_data)[nca].next = prev;
+
+ for (int i = k - 2; i >= 0; i -= 2) {
+ subblossoms.push_back(right_path[i + 1]);
+ (*_blossom_data)[right_path[i + 1]].status = EVEN;
+ oddToEven(right_path[i + 1], tree);
+ _tree_set->erase(right_path[i + 1]);
+
+ (*_blossom_data)[right_path[i + 1]].next =
+ (*_blossom_data)[right_path[i + 1]].pred;
+
+ subblossoms.push_back(right_path[i]);
+ _tree_set->erase(right_path[i]);
+ }
+
+ int surface =
+ _blossom_set->join(subblossoms.begin(), subblossoms.end());
+
+ for (int i = 0; i < int(subblossoms.size()); ++i) {
+ if (!_blossom_set->trivial(subblossoms[i])) {
+ (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+ }
+ (*_blossom_data)[subblossoms[i]].status = MATCHED;
+ }
+
+ (*_blossom_data)[surface].pot = -2 * _delta_sum;
+ (*_blossom_data)[surface].offset = 0;
+ (*_blossom_data)[surface].status = EVEN;
+ (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
+ (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
+
+ _tree_set->insert(surface, tree);
+ _tree_set->erase(nca);
+ }
+
+ void splitBlossom(int blossom) {
+ Arc next = (*_blossom_data)[blossom].next;
+ Arc pred = (*_blossom_data)[blossom].pred;
+
+ int tree = _tree_set->find(blossom);
+
+ (*_blossom_data)[blossom].status = MATCHED;
+ oddToMatched(blossom);
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+
+ std::vector subblossoms;
+ _blossom_set->split(blossom, std::back_inserter(subblossoms));
+
+ Value offset = (*_blossom_data)[blossom].offset;
+ int b = _blossom_set->find(_graph.source(pred));
+ int d = _blossom_set->find(_graph.source(next));
+
+ int ib = -1, id = -1;
+ for (int i = 0; i < int(subblossoms.size()); ++i) {
+ if (subblossoms[i] == b) ib = i;
+ if (subblossoms[i] == d) id = i;
+
+ (*_blossom_data)[subblossoms[i]].offset = offset;
+ if (!_blossom_set->trivial(subblossoms[i])) {
+ (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+ }
+ if (_blossom_set->classPrio(subblossoms[i]) !=
+ std::numeric_limits::max()) {
+ _delta2->push(subblossoms[i],
+ _blossom_set->classPrio(subblossoms[i]) -
+ (*_blossom_data)[subblossoms[i]].offset);
+ }
+ }
+
+ if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
+ for (int i = (id + 1) % subblossoms.size();
+ i != ib; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ (*_blossom_data)[sb].next =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+ }
+
+ for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+ (*_blossom_data)[sb].status = ODD;
+ matchedToOdd(sb);
+ _tree_set->insert(sb, tree);
+ (*_blossom_data)[sb].pred = pred;
+ (*_blossom_data)[sb].next =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+
+ pred = (*_blossom_data)[ub].next;
+
+ (*_blossom_data)[tb].status = EVEN;
+ matchedToEven(tb, tree);
+ _tree_set->insert(tb, tree);
+ (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+ }
+
+ (*_blossom_data)[subblossoms[id]].status = ODD;
+ matchedToOdd(subblossoms[id]);
+ _tree_set->insert(subblossoms[id], tree);
+ (*_blossom_data)[subblossoms[id]].next = next;
+ (*_blossom_data)[subblossoms[id]].pred = pred;
+
+ } else {
+
+ for (int i = (ib + 1) % subblossoms.size();
+ i != id; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ (*_blossom_data)[sb].next =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+ }
+
+ for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+ (*_blossom_data)[sb].status = ODD;
+ matchedToOdd(sb);
+ _tree_set->insert(sb, tree);
+ (*_blossom_data)[sb].next = next;
+ (*_blossom_data)[sb].pred =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+
+ (*_blossom_data)[tb].status = EVEN;
+ matchedToEven(tb, tree);
+ _tree_set->insert(tb, tree);
+ (*_blossom_data)[tb].pred =
+ (*_blossom_data)[tb].next =
+ _graph.oppositeArc((*_blossom_data)[ub].next);
+ next = (*_blossom_data)[ub].next;
+ }
+
+ (*_blossom_data)[subblossoms[ib]].status = ODD;
+ matchedToOdd(subblossoms[ib]);
+ _tree_set->insert(subblossoms[ib], tree);
+ (*_blossom_data)[subblossoms[ib]].next = next;
+ (*_blossom_data)[subblossoms[ib]].pred = pred;
+ }
+ _tree_set->erase(blossom);
+ }
+
+ void extractBlossom(int blossom, const Node& base, const Arc& matching) {
+ if (_blossom_set->trivial(blossom)) {
+ int bi = (*_node_index)[base];
+ Value pot = (*_node_data)[bi].pot;
+
+ _matching->set(base, matching);
+ _blossom_node_list.push_back(base);
+ _node_potential->set(base, pot);
+ } else {
+
+ Value pot = (*_blossom_data)[blossom].pot;
+ int bn = _blossom_node_list.size();
+
+ std::vector subblossoms;
+ _blossom_set->split(blossom, std::back_inserter(subblossoms));
+ int b = _blossom_set->find(base);
+ int ib = -1;
+ for (int i = 0; i < int(subblossoms.size()); ++i) {
+ if (subblossoms[i] == b) { ib = i; break; }
+ }
+
+ for (int i = 1; i < int(subblossoms.size()); i += 2) {
+ int sb = subblossoms[(ib + i) % subblossoms.size()];
+ int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
+
+ Arc m = (*_blossom_data)[tb].next;
+ extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
+ extractBlossom(tb, _graph.source(m), m);
+ }
+ extractBlossom(subblossoms[ib], base, matching);
+
+ int en = _blossom_node_list.size();
+
+ _blossom_potential.push_back(BlossomVariable(bn, en, pot));
+ }
+ }
+
+ void extractMatching() {
+ std::vector blossoms;
+ for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
+ blossoms.push_back(c);
+ }
+
+ for (int i = 0; i < int(blossoms.size()); ++i) {
+ if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
+
+ Value offset = (*_blossom_data)[blossoms[i]].offset;
+ (*_blossom_data)[blossoms[i]].pot += 2 * offset;
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
+ n != INVALID; ++n) {
+ (*_node_data)[(*_node_index)[n]].pot -= offset;
+ }
+
+ Arc matching = (*_blossom_data)[blossoms[i]].next;
+ Node base = _graph.source(matching);
+ extractBlossom(blossoms[i], base, matching);
+ } else {
+ Node base = (*_blossom_data)[blossoms[i]].base;
+ extractBlossom(blossoms[i], base, INVALID);
+ }
+ }
+ }
+
+ public:
+
+ /// \brief Constructor
+ ///
+ /// Constructor.
+ MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
+ : _graph(graph), _weight(weight), _matching(0),
+ _node_potential(0), _blossom_potential(), _blossom_node_list(),
+ _node_num(0), _blossom_num(0),
+
+ _blossom_index(0), _blossom_set(0), _blossom_data(0),
+ _node_index(0), _node_heap_index(0), _node_data(0),
+ _tree_set_index(0), _tree_set(0),
+
+ _delta1_index(0), _delta1(0),
+ _delta2_index(0), _delta2(0),
+ _delta3_index(0), _delta3(0),
+ _delta4_index(0), _delta4(0),
+
+ _delta_sum() {}
+
+ ~MaxWeightedMatching() {
+ destroyStructures();
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to use the
+ /// \c run() member function.
+
+ ///@{
+
+ /// \brief Initialize the algorithm
+ ///
+ /// Initialize the algorithm
+ void init() {
+ createStructures();
+
+ for (ArcIt e(_graph); e != INVALID; ++e) {
+ _node_heap_index->set(e, BinHeap::PRE_HEAP);
+ }
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _delta1_index->set(n, _delta1->PRE_HEAP);
+ }
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _delta3_index->set(e, _delta3->PRE_HEAP);
+ }
+ for (int i = 0; i < _blossom_num; ++i) {
+ _delta2_index->set(i, _delta2->PRE_HEAP);
+ _delta4_index->set(i, _delta4->PRE_HEAP);
+ }
+
+ int index = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ Value max = 0;
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ if (_graph.target(e) == n) continue;
+ if ((dualScale * _weight[e]) / 2 > max) {
+ max = (dualScale * _weight[e]) / 2;
+ }
+ }
+ _node_index->set(n, index);
+ (*_node_data)[index].pot = max;
+ _delta1->push(n, max);
+ int blossom =
+ _blossom_set->insert(n, std::numeric_limits::max());
+
+ _tree_set->insert(blossom);
+
+ (*_blossom_data)[blossom].status = EVEN;
+ (*_blossom_data)[blossom].pred = INVALID;
+ (*_blossom_data)[blossom].next = INVALID;
+ (*_blossom_data)[blossom].pot = 0;
+ (*_blossom_data)[blossom].offset = 0;
+ ++index;
+ }
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ int si = (*_node_index)[_graph.u(e)];
+ int ti = (*_node_index)[_graph.v(e)];
+ if (_graph.u(e) != _graph.v(e)) {
+ _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
+ dualScale * _weight[e]) / 2);
+ }
+ }
+ }
+
+ /// \brief Starts the algorithm
+ ///
+ /// Starts the algorithm
+ void start() {
+ enum OpType {
+ D1, D2, D3, D4
+ };
+
+ int unmatched = _node_num;
+ while (unmatched > 0) {
+ Value d1 = !_delta1->empty() ?
+ _delta1->prio() : std::numeric_limits::max();
+
+ Value d2 = !_delta2->empty() ?
+ _delta2->prio() : std::numeric_limits::max();
+
+ Value d3 = !_delta3->empty() ?
+ _delta3->prio() : std::numeric_limits::max();
+
+ Value d4 = !_delta4->empty() ?
+ _delta4->prio() : std::numeric_limits::max();
+
+ _delta_sum = d1; OpType ot = D1;
+ if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
+ if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
+ if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
+
+
+ switch (ot) {
+ case D1:
+ {
+ Node n = _delta1->top();
+ unmatchNode(n);
+ --unmatched;
+ }
+ break;
+ case D2:
+ {
+ int blossom = _delta2->top();
+ Node n = _blossom_set->classTop(blossom);
+ Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
+ extendOnArc(e);
+ }
+ break;
+ case D3:
+ {
+ Edge e = _delta3->top();
+
+ int left_blossom = _blossom_set->find(_graph.u(e));
+ int right_blossom = _blossom_set->find(_graph.v(e));
+
+ if (left_blossom == right_blossom) {
+ _delta3->pop();
+ } else {
+ int left_tree;
+ if ((*_blossom_data)[left_blossom].status == EVEN) {
+ left_tree = _tree_set->find(left_blossom);
+ } else {
+ left_tree = -1;
+ ++unmatched;
+ }
+ int right_tree;
+ if ((*_blossom_data)[right_blossom].status == EVEN) {
+ right_tree = _tree_set->find(right_blossom);
+ } else {
+ right_tree = -1;
+ ++unmatched;
+ }
+
+ if (left_tree == right_tree) {
+ shrinkOnEdge(e, left_tree);
+ } else {
+ augmentOnEdge(e);
+ unmatched -= 2;
+ }
+ }
+ } break;
+ case D4:
+ splitBlossom(_delta4->top());
+ break;
+ }
+ }
+ extractMatching();
+ }
+
+ /// \brief Runs %MaxWeightedMatching algorithm.
+ ///
+ /// This method runs the %MaxWeightedMatching algorithm.
+ ///
+ /// \note mwm.run() is just a shortcut of the following code.
+ /// \code
+ /// mwm.init();
+ /// mwm.start();
+ /// \endcode
+ void run() {
+ init();
+ start();
+ }
+
+ /// @}
+
+ /// \name Primal solution
+ /// Functions to get the primal solution, ie. the matching.
+
+ /// @{
+
+ /// \brief Returns the weight of the matching.
+ ///
+ /// Returns the weight of the matching.
+ Value matchingValue() const {
+ Value sum = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_matching)[n] != INVALID) {
+ sum += _weight[(*_matching)[n]];
+ }
+ }
+ return sum /= 2;
+ }
+
+ /// \brief Returns the cardinality of the matching.
+ ///
+ /// Returns the cardinality of the matching.
+ int matchingSize() const {
+ int num = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_matching)[n] != INVALID) {
+ ++num;
+ }
+ }
+ return num /= 2;
+ }
+
+ /// \brief Returns true when the edge is in the matching.
+ ///
+ /// Returns true when the edge is in the matching.
+ bool matching(const Edge& edge) const {
+ return edge == (*_matching)[_graph.u(edge)];
+ }
+
+ /// \brief Returns the incident matching arc.
+ ///
+ /// Returns the incident matching arc from given node. If the
+ /// node is not matched then it gives back \c INVALID.
+ Arc matching(const Node& node) const {
+ return (*_matching)[node];
+ }
+
+ /// \brief Returns the mate of the node.
+ ///
+ /// Returns the adjancent node in a mathcing arc. If the node is
+ /// not matched then it gives back \c INVALID.
+ Node mate(const Node& node) const {
+ return (*_matching)[node] != INVALID ?
+ _graph.target((*_matching)[node]) : INVALID;
+ }
+
+ /// @}
+
+ /// \name Dual solution
+ /// Functions to get the dual solution.
+
+ /// @{
+
+ /// \brief Returns the value of the dual solution.
+ ///
+ /// Returns the value of the dual solution. It should be equal to
+ /// the primal value scaled by \ref dualScale "dual scale".
+ Value dualValue() const {
+ Value sum = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ sum += nodeValue(n);
+ }
+ for (int i = 0; i < blossomNum(); ++i) {
+ sum += blossomValue(i) * (blossomSize(i) / 2);
+ }
+ return sum;
+ }
+
+ /// \brief Returns the value of the node.
+ ///
+ /// Returns the the value of the node.
+ Value nodeValue(const Node& n) const {
+ return (*_node_potential)[n];
+ }
+
+ /// \brief Returns the number of the blossoms in the basis.
+ ///
+ /// Returns the number of the blossoms in the basis.
+ /// \see BlossomIt
+ int blossomNum() const {
+ return _blossom_potential.size();
+ }
+
+
+ /// \brief Returns the number of the nodes in the blossom.
+ ///
+ /// Returns the number of the nodes in the blossom.
+ int blossomSize(int k) const {
+ return _blossom_potential[k].end - _blossom_potential[k].begin;
+ }
+
+ /// \brief Returns the value of the blossom.
+ ///
+ /// Returns the the value of the blossom.
+ /// \see BlossomIt
+ Value blossomValue(int k) const {
+ return _blossom_potential[k].value;
+ }
+
+ /// \brief Iterator for obtaining the nodes of the blossom.
+ ///
+ /// Iterator for obtaining the nodes of the blossom. This class
+ /// provides a common lemon style iterator for listing a
+ /// subset of the nodes.
+ class BlossomIt {
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// Constructor to get the nodes of the variable.
+ BlossomIt(const MaxWeightedMatching& algorithm, int variable)
+ : _algorithm(&algorithm)
+ {
+ _index = _algorithm->_blossom_potential[variable].begin;
+ _last = _algorithm->_blossom_potential[variable].end;
+ }
+
+ /// \brief Conversion to node.
+ ///
+ /// Conversion to node.
+ operator Node() const {
+ return _algorithm->_blossom_node_list[_index];
+ }
+
+ /// \brief Increment operator.
+ ///
+ /// Increment operator.
+ BlossomIt& operator++() {
+ ++_index;
+ return *this;
+ }
+
+ /// \brief Validity checking
+ ///
+ /// Checks whether the iterator is invalid.
+ bool operator==(Invalid) const { return _index == _last; }
+
+ /// \brief Validity checking
+ ///
+ /// Checks whether the iterator is valid.
+ bool operator!=(Invalid) const { return _index != _last; }
+
+ private:
+ const MaxWeightedMatching* _algorithm;
+ int _last;
+ int _index;
+ };
+
+ /// @}
+
+ };
+
+ /// \ingroup matching
+ ///
+ /// \brief Weighted perfect matching in general graphs
+ ///
+ /// This class provides an efficient implementation of Edmond's
+ /// maximum weighted perfect matching algorithm. The implementation
+ /// is based on extensive use of priority queues and provides
+ /// \f$O(nm\log(n))\f$ time complexity.
+ ///
+ /// The maximum weighted matching problem is to find undirected
+ /// edges in the graph with maximum overall weight and no two of
+ /// them shares their ends and covers all nodes. The problem can be
+ /// formulated with the following linear program.
+ /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
+ /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
+ \quad \forall B\in\mathcal{O}\f] */
+ /// \f[x_e \ge 0\quad \forall e\in E\f]
+ /// \f[\max \sum_{e\in E}x_ew_e\f]
+ /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+ /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
+ /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
+ /// subsets of the nodes.
+ ///
+ /// The algorithm calculates an optimal matching and a proof of the
+ /// optimality. The solution of the dual problem can be used to check
+ /// the result of the algorithm. The dual linear problem is the
+ /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
+ w_{uv} \quad \forall uv\in E\f] */
+ /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
+ /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
+ \frac{\vert B \vert - 1}{2}z_B\f] */
+ ///
+ /// The algorithm can be executed with \c run() or the \c init() and
+ /// then the \c start() member functions. After it the matching can
+ /// be asked with \c matching() or mate() functions. The dual
+ /// solution can be get with \c nodeValue(), \c blossomNum() and \c
+ /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
+ /// "BlossomIt" nested class which is able to iterate on the nodes
+ /// of a blossom. If the value type is integral then the dual
+ /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
+ template >
+ class MaxWeightedPerfectMatching {
+ public:
+
+ typedef _Graph Graph;
+ typedef _WeightMap WeightMap;
+ typedef typename WeightMap::Value Value;
+
+ /// \brief Scaling factor for dual solution
+ ///
+ /// Scaling factor for dual solution, it is equal to 4 or 1
+ /// according to the value type.
+ static const int dualScale =
+ std::numeric_limits::is_integer ? 4 : 1;
+
+ typedef typename Graph::template NodeMap
+ MatchingMap;
+
+ private:
+
+ TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+ typedef typename Graph::template NodeMap NodePotential;
+ typedef std::vector BlossomNodeList;
+
+ struct BlossomVariable {
+ int begin, end;
+ Value value;
+
+ BlossomVariable(int _begin, int _end, Value _value)
+ : begin(_begin), end(_end), value(_value) {}
+
+ };
+
+ typedef std::vector BlossomPotential;
+
+ const Graph& _graph;
+ const WeightMap& _weight;
+
+ MatchingMap* _matching;
+
+ NodePotential* _node_potential;
+
+ BlossomPotential _blossom_potential;
+ BlossomNodeList _blossom_node_list;
+
+ int _node_num;
+ int _blossom_num;
+
+ typedef RangeMap IntIntMap;
+
+ enum Status {
+ EVEN = -1, MATCHED = 0, ODD = 1
+ };
+
+ typedef HeapUnionFind BlossomSet;
+ struct BlossomData {
+ int tree;
+ Status status;
+ Arc pred, next;
+ Value pot, offset;
+ };
+
+ IntNodeMap *_blossom_index;
+ BlossomSet *_blossom_set;
+ RangeMap* _blossom_data;
+
+ IntNodeMap *_node_index;
+ IntArcMap *_node_heap_index;
+
+ struct NodeData {
+
+ NodeData(IntArcMap& node_heap_index)
+ : heap(node_heap_index) {}
+
+ int blossom;
+ Value pot;
+ BinHeap heap;
+ std::map heap_index;
+
+ int tree;
+ };
+
+ RangeMap* _node_data;
+
+ typedef ExtendFindEnum TreeSet;
+
+ IntIntMap *_tree_set_index;
+ TreeSet *_tree_set;
+
+ IntIntMap *_delta2_index;
+ BinHeap *_delta2;
+
+ IntEdgeMap *_delta3_index;
+ BinHeap *_delta3;
+
+ IntIntMap *_delta4_index;
+ BinHeap *_delta4;
+
+ Value _delta_sum;
+
+ void createStructures() {
+ _node_num = countNodes(_graph);
+ _blossom_num = _node_num * 3 / 2;
+
+ if (!_matching) {
+ _matching = new MatchingMap(_graph);
+ }
+ if (!_node_potential) {
+ _node_potential = new NodePotential(_graph);
+ }
+ if (!_blossom_set) {
+ _blossom_index = new IntNodeMap(_graph);
+ _blossom_set = new BlossomSet(*_blossom_index);
+ _blossom_data = new RangeMap(_blossom_num);
+ }
+
+ if (!_node_index) {
+ _node_index = new IntNodeMap(_graph);
+ _node_heap_index = new IntArcMap(_graph);
+ _node_data = new RangeMap(_node_num,
+ NodeData(*_node_heap_index));
+ }
+
+ if (!_tree_set) {
+ _tree_set_index = new IntIntMap(_blossom_num);
+ _tree_set = new TreeSet(*_tree_set_index);
+ }
+ if (!_delta2) {
+ _delta2_index = new IntIntMap(_blossom_num);
+ _delta2 = new BinHeap(*_delta2_index);
+ }
+ if (!_delta3) {
+ _delta3_index = new IntEdgeMap(_graph);
+ _delta3 = new BinHeap(*_delta3_index);
+ }
+ if (!_delta4) {
+ _delta4_index = new IntIntMap(_blossom_num);
+ _delta4 = new BinHeap(*_delta4_index);
+ }
+ }
+
+ void destroyStructures() {
+ _node_num = countNodes(_graph);
+ _blossom_num = _node_num * 3 / 2;
+
+ if (_matching) {
+ delete _matching;
+ }
+ if (_node_potential) {
+ delete _node_potential;
+ }
+ if (_blossom_set) {
+ delete _blossom_index;
+ delete _blossom_set;
+ delete _blossom_data;
+ }
+
+ if (_node_index) {
+ delete _node_index;
+ delete _node_heap_index;
+ delete _node_data;
+ }
+
+ if (_tree_set) {
+ delete _tree_set_index;
+ delete _tree_set;
+ }
+ if (_delta2) {
+ delete _delta2_index;
+ delete _delta2;
+ }
+ if (_delta3) {
+ delete _delta3_index;
+ delete _delta3;
+ }
+ if (_delta4) {
+ delete _delta4_index;
+ delete _delta4;
+ }
+ }
+
+ void matchedToEven(int blossom, int tree) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+
+ if (!_blossom_set->trivial(blossom)) {
+ (*_blossom_data)[blossom].pot -=
+ 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+
+ _blossom_set->increase(n, std::numeric_limits::max());
+ int ni = (*_node_index)[n];
+
+ (*_node_data)[ni].heap.clear();
+ (*_node_data)[ni].heap_index.clear();
+
+ (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if ((*_blossom_data)[vb].status == EVEN) {
+ if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+ _delta3->push(e, rw / 2);
+ }
+ } else {
+ typename std::map::iterator it =
+ (*_node_data)[vi].heap_index.find(tree);
+
+ if (it != (*_node_data)[vi].heap_index.end()) {
+ if ((*_node_data)[vi].heap[it->second] > rw) {
+ (*_node_data)[vi].heap.replace(it->second, e);
+ (*_node_data)[vi].heap.decrease(e, rw);
+ it->second = e;
+ }
+ } else {
+ (*_node_data)[vi].heap.push(e, rw);
+ (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+ }
+
+ if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+ _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+ if ((*_blossom_data)[vb].status == MATCHED) {
+ if (_delta2->state(vb) != _delta2->IN_HEAP) {
+ _delta2->push(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset){
+ _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ }
+ }
+ }
+ }
+ }
+ }
+ (*_blossom_data)[blossom].offset = 0;
+ }
+
+ void matchedToOdd(int blossom) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+ (*_blossom_data)[blossom].offset += _delta_sum;
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+
+ void evenToMatched(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ (*_blossom_data)[blossom].pot += 2 * _delta_sum;
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+ (*_node_data)[ni].pot -= _delta_sum;
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if (vb == blossom) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else if ((*_blossom_data)[vb].status == EVEN) {
+
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+
+ int vt = _tree_set->find(vb);
+
+ if (vt != tree) {
+
+ Arc r = _graph.oppositeArc(e);
+
+ typename std::map::iterator it =
+ (*_node_data)[ni].heap_index.find(vt);
+
+ if (it != (*_node_data)[ni].heap_index.end()) {
+ if ((*_node_data)[ni].heap[it->second] > rw) {
+ (*_node_data)[ni].heap.replace(it->second, r);
+ (*_node_data)[ni].heap.decrease(r, rw);
+ it->second = r;
+ }
+ } else {
+ (*_node_data)[ni].heap.push(r, rw);
+ (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+ }
+
+ if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
+ _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+
+ if (_delta2->state(blossom) != _delta2->IN_HEAP) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ } else if ((*_delta2)[blossom] >
+ _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset){
+ _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+ }
+ } else {
+
+ typename std::map::iterator it =
+ (*_node_data)[vi].heap_index.find(tree);
+
+ if (it != (*_node_data)[vi].heap_index.end()) {
+ (*_node_data)[vi].heap.erase(it->second);
+ (*_node_data)[vi].heap_index.erase(it);
+ if ((*_node_data)[vi].heap.empty()) {
+ _blossom_set->increase(v, std::numeric_limits::max());
+ } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
+ _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
+ }
+
+ if ((*_blossom_data)[vb].status == MATCHED) {
+ if (_blossom_set->classPrio(vb) ==
+ std::numeric_limits::max()) {
+ _delta2->erase(vb);
+ } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset) {
+ _delta2->increase(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void oddToMatched(int blossom) {
+ (*_blossom_data)[blossom].offset -= _delta_sum;
+
+ if (_blossom_set->classPrio(blossom) !=
+ std::numeric_limits::max()) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ }
+ }
+
+ void oddToEven(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ (*_blossom_data)[blossom].pot -=
+ 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+
+ _blossom_set->increase(n, std::numeric_limits::max());
+
+ (*_node_data)[ni].heap.clear();
+ (*_node_data)[ni].heap_index.clear();
+ (*_node_data)[ni].pot +=
+ 2 * _delta_sum - (*_blossom_data)[blossom].offset;
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.source(e);
+ int vb = _blossom_set->find(v);
+ int vi = (*_node_index)[v];
+
+ Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+ dualScale * _weight[e];
+
+ if ((*_blossom_data)[vb].status == EVEN) {
+ if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+ _delta3->push(e, rw / 2);
+ }
+ } else {
+
+ typename std::map::iterator it =
+ (*_node_data)[vi].heap_index.find(tree);
+
+ if (it != (*_node_data)[vi].heap_index.end()) {
+ if ((*_node_data)[vi].heap[it->second] > rw) {
+ (*_node_data)[vi].heap.replace(it->second, e);
+ (*_node_data)[vi].heap.decrease(e, rw);
+ it->second = e;
+ }
+ } else {
+ (*_node_data)[vi].heap.push(e, rw);
+ (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+ }
+
+ if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+ _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+ if ((*_blossom_data)[vb].status == MATCHED) {
+ if (_delta2->state(vb) != _delta2->IN_HEAP) {
+ _delta2->push(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset) {
+ _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+ (*_blossom_data)[vb].offset);
+ }
+ }
+ }
+ }
+ }
+ }
+ (*_blossom_data)[blossom].offset = 0;
+ }
+
+ void alternatePath(int even, int tree) {
+ int odd;
+
+ evenToMatched(even, tree);
+ (*_blossom_data)[even].status = MATCHED;
+
+ while ((*_blossom_data)[even].pred != INVALID) {
+ odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
+ (*_blossom_data)[odd].status = MATCHED;
+ oddToMatched(odd);
+ (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
+
+ even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
+ (*_blossom_data)[even].status = MATCHED;
+ evenToMatched(even, tree);
+ (*_blossom_data)[even].next =
+ _graph.oppositeArc((*_blossom_data)[odd].pred);
+ }
+
+ }
+
+ void destroyTree(int tree) {
+ for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
+ if ((*_blossom_data)[b].status == EVEN) {
+ (*_blossom_data)[b].status = MATCHED;
+ evenToMatched(b, tree);
+ } else if ((*_blossom_data)[b].status == ODD) {
+ (*_blossom_data)[b].status = MATCHED;
+ oddToMatched(b);
+ }
+ }
+ _tree_set->eraseClass(tree);
+ }
+
+ void augmentOnEdge(const Edge& edge) {
+
+ int left = _blossom_set->find(_graph.u(edge));
+ int right = _blossom_set->find(_graph.v(edge));
+
+ int left_tree = _tree_set->find(left);
+ alternatePath(left, left_tree);
+ destroyTree(left_tree);
+
+ int right_tree = _tree_set->find(right);
+ alternatePath(right, right_tree);
+ destroyTree(right_tree);
+
+ (*_blossom_data)[left].next = _graph.direct(edge, true);
+ (*_blossom_data)[right].next = _graph.direct(edge, false);
+ }
+
+ void extendOnArc(const Arc& arc) {
+ int base = _blossom_set->find(_graph.target(arc));
+ int tree = _tree_set->find(base);
+
+ int odd = _blossom_set->find(_graph.source(arc));
+ _tree_set->insert(odd, tree);
+ (*_blossom_data)[odd].status = ODD;
+ matchedToOdd(odd);
+ (*_blossom_data)[odd].pred = arc;
+
+ int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
+ (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
+ _tree_set->insert(even, tree);
+ (*_blossom_data)[even].status = EVEN;
+ matchedToEven(even, tree);
+ }
+
+ void shrinkOnEdge(const Edge& edge, int tree) {
+ int nca = -1;
+ std::vector left_path, right_path;
+
+ {
+ std::set left_set, right_set;
+ int left = _blossom_set->find(_graph.u(edge));
+ left_path.push_back(left);
+ left_set.insert(left);
+
+ int right = _blossom_set->find(_graph.v(edge));
+ right_path.push_back(right);
+ right_set.insert(right);
+
+ while (true) {
+
+ if ((*_blossom_data)[left].pred == INVALID) break;
+
+ left =
+ _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
+ left_path.push_back(left);
+ left =
+ _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
+ left_path.push_back(left);
+
+ left_set.insert(left);
+
+ if (right_set.find(left) != right_set.end()) {
+ nca = left;
+ break;
+ }
+
+ if ((*_blossom_data)[right].pred == INVALID) break;
+
+ right =
+ _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
+ right_path.push_back(right);
+ right =
+ _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
+ right_path.push_back(right);
+
+ right_set.insert(right);
+
+ if (left_set.find(right) != left_set.end()) {
+ nca = right;
+ break;
+ }
+
+ }
+
+ if (nca == -1) {
+ if ((*_blossom_data)[left].pred == INVALID) {
+ nca = right;
+ while (left_set.find(nca) == left_set.end()) {
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ right_path.push_back(nca);
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ right_path.push_back(nca);
+ }
+ } else {
+ nca = left;
+ while (right_set.find(nca) == right_set.end()) {
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ nca =
+ _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ }
+ }
+ }
+ }
+
+ std::vector subblossoms;
+ Arc prev;
+
+ prev = _graph.direct(edge, true);
+ for (int i = 0; left_path[i] != nca; i += 2) {
+ subblossoms.push_back(left_path[i]);
+ (*_blossom_data)[left_path[i]].next = prev;
+ _tree_set->erase(left_path[i]);
+
+ subblossoms.push_back(left_path[i + 1]);
+ (*_blossom_data)[left_path[i + 1]].status = EVEN;
+ oddToEven(left_path[i + 1], tree);
+ _tree_set->erase(left_path[i + 1]);
+ prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
+ }
+
+ int k = 0;
+ while (right_path[k] != nca) ++k;
+
+ subblossoms.push_back(nca);
+ (*_blossom_data)[nca].next = prev;
+
+ for (int i = k - 2; i >= 0; i -= 2) {
+ subblossoms.push_back(right_path[i + 1]);
+ (*_blossom_data)[right_path[i + 1]].status = EVEN;
+ oddToEven(right_path[i + 1], tree);
+ _tree_set->erase(right_path[i + 1]);
+
+ (*_blossom_data)[right_path[i + 1]].next =
+ (*_blossom_data)[right_path[i + 1]].pred;
+
+ subblossoms.push_back(right_path[i]);
+ _tree_set->erase(right_path[i]);
+ }
+
+ int surface =
+ _blossom_set->join(subblossoms.begin(), subblossoms.end());
+
+ for (int i = 0; i < int(subblossoms.size()); ++i) {
+ if (!_blossom_set->trivial(subblossoms[i])) {
+ (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+ }
+ (*_blossom_data)[subblossoms[i]].status = MATCHED;
+ }
+
+ (*_blossom_data)[surface].pot = -2 * _delta_sum;
+ (*_blossom_data)[surface].offset = 0;
+ (*_blossom_data)[surface].status = EVEN;
+ (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
+ (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
+
+ _tree_set->insert(surface, tree);
+ _tree_set->erase(nca);
+ }
+
+ void splitBlossom(int blossom) {
+ Arc next = (*_blossom_data)[blossom].next;
+ Arc pred = (*_blossom_data)[blossom].pred;
+
+ int tree = _tree_set->find(blossom);
+
+ (*_blossom_data)[blossom].status = MATCHED;
+ oddToMatched(blossom);
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+
+ std::vector subblossoms;
+ _blossom_set->split(blossom, std::back_inserter(subblossoms));
+
+ Value offset = (*_blossom_data)[blossom].offset;
+ int b = _blossom_set->find(_graph.source(pred));
+ int d = _blossom_set->find(_graph.source(next));
+
+ int ib = -1, id = -1;
+ for (int i = 0; i < int(subblossoms.size()); ++i) {
+ if (subblossoms[i] == b) ib = i;
+ if (subblossoms[i] == d) id = i;
+
+ (*_blossom_data)[subblossoms[i]].offset = offset;
+ if (!_blossom_set->trivial(subblossoms[i])) {
+ (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+ }
+ if (_blossom_set->classPrio(subblossoms[i]) !=
+ std::numeric_limits::max()) {
+ _delta2->push(subblossoms[i],
+ _blossom_set->classPrio(subblossoms[i]) -
+ (*_blossom_data)[subblossoms[i]].offset);
+ }
+ }
+
+ if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
+ for (int i = (id + 1) % subblossoms.size();
+ i != ib; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ (*_blossom_data)[sb].next =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+ }
+
+ for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+ (*_blossom_data)[sb].status = ODD;
+ matchedToOdd(sb);
+ _tree_set->insert(sb, tree);
+ (*_blossom_data)[sb].pred = pred;
+ (*_blossom_data)[sb].next =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+
+ pred = (*_blossom_data)[ub].next;
+
+ (*_blossom_data)[tb].status = EVEN;
+ matchedToEven(tb, tree);
+ _tree_set->insert(tb, tree);
+ (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+ }
+
+ (*_blossom_data)[subblossoms[id]].status = ODD;
+ matchedToOdd(subblossoms[id]);
+ _tree_set->insert(subblossoms[id], tree);
+ (*_blossom_data)[subblossoms[id]].next = next;
+ (*_blossom_data)[subblossoms[id]].pred = pred;
+
+ } else {
+
+ for (int i = (ib + 1) % subblossoms.size();
+ i != id; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ (*_blossom_data)[sb].next =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+ }
+
+ for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
+ int sb = subblossoms[i];
+ int tb = subblossoms[(i + 1) % subblossoms.size()];
+ int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+ (*_blossom_data)[sb].status = ODD;
+ matchedToOdd(sb);
+ _tree_set->insert(sb, tree);
+ (*_blossom_data)[sb].next = next;
+ (*_blossom_data)[sb].pred =
+ _graph.oppositeArc((*_blossom_data)[tb].next);
+
+ (*_blossom_data)[tb].status = EVEN;
+ matchedToEven(tb, tree);
+ _tree_set->insert(tb, tree);
+ (*_blossom_data)[tb].pred =
+ (*_blossom_data)[tb].next =
+ _graph.oppositeArc((*_blossom_data)[ub].next);
+ next = (*_blossom_data)[ub].next;
+ }
+
+ (*_blossom_data)[subblossoms[ib]].status = ODD;
+ matchedToOdd(subblossoms[ib]);
+ _tree_set->insert(subblossoms[ib], tree);
+ (*_blossom_data)[subblossoms[ib]].next = next;
+ (*_blossom_data)[subblossoms[ib]].pred = pred;
+ }
+ _tree_set->erase(blossom);
+ }
+
+ void extractBlossom(int blossom, const Node& base, const Arc& matching) {
+ if (_blossom_set->trivial(blossom)) {
+ int bi = (*_node_index)[base];
+ Value pot = (*_node_data)[bi].pot;
+
+ _matching->set(base, matching);
+ _blossom_node_list.push_back(base);
+ _node_potential->set(base, pot);
+ } else {
+
+ Value pot = (*_blossom_data)[blossom].pot;
+ int bn = _blossom_node_list.size();
+
+ std::vector subblossoms;
+ _blossom_set->split(blossom, std::back_inserter(subblossoms));
+ int b = _blossom_set->find(base);
+ int ib = -1;
+ for (int i = 0; i < int(subblossoms.size()); ++i) {
+ if (subblossoms[i] == b) { ib = i; break; }
+ }
+
+ for (int i = 1; i < int(subblossoms.size()); i += 2) {
+ int sb = subblossoms[(ib + i) % subblossoms.size()];
+ int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
+
+ Arc m = (*_blossom_data)[tb].next;
+ extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
+ extractBlossom(tb, _graph.source(m), m);
+ }
+ extractBlossom(subblossoms[ib], base, matching);
+
+ int en = _blossom_node_list.size();
+
+ _blossom_potential.push_back(BlossomVariable(bn, en, pot));
+ }
+ }
+
+ void extractMatching() {
+ std::vector blossoms;
+ for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
+ blossoms.push_back(c);
+ }
+
+ for (int i = 0; i < int(blossoms.size()); ++i) {
+
+ Value offset = (*_blossom_data)[blossoms[i]].offset;
+ (*_blossom_data)[blossoms[i]].pot += 2 * offset;
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
+ n != INVALID; ++n) {
+ (*_node_data)[(*_node_index)[n]].pot -= offset;
+ }
+
+ Arc matching = (*_blossom_data)[blossoms[i]].next;
+ Node base = _graph.source(matching);
+ extractBlossom(blossoms[i], base, matching);
+ }
+ }
+
+ public:
+
+ /// \brief Constructor
+ ///
+ /// Constructor.
+ MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
+ : _graph(graph), _weight(weight), _matching(0),
+ _node_potential(0), _blossom_potential(), _blossom_node_list(),
+ _node_num(0), _blossom_num(0),
+
+ _blossom_index(0), _blossom_set(0), _blossom_data(0),
+ _node_index(0), _node_heap_index(0), _node_data(0),
+ _tree_set_index(0), _tree_set(0),
+
+ _delta2_index(0), _delta2(0),
+ _delta3_index(0), _delta3(0),
+ _delta4_index(0), _delta4(0),
+
+ _delta_sum() {}
+
+ ~MaxWeightedPerfectMatching() {
+ destroyStructures();
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to use the
+ /// \c run() member function.
+
+ ///@{
+
+ /// \brief Initialize the algorithm
+ ///
+ /// Initialize the algorithm
+ void init() {
+ createStructures();
+
+ for (ArcIt e(_graph); e != INVALID; ++e) {
+ _node_heap_index->set(e, BinHeap::PRE_HEAP);
+ }
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ _delta3_index->set(e, _delta3->PRE_HEAP);
+ }
+ for (int i = 0; i < _blossom_num; ++i) {
+ _delta2_index->set(i, _delta2->PRE_HEAP);
+ _delta4_index->set(i, _delta4->PRE_HEAP);
+ }
+
+ int index = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ Value max = - std::numeric_limits::max();
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ if (_graph.target(e) == n) continue;
+ if ((dualScale * _weight[e]) / 2 > max) {
+ max = (dualScale * _weight[e]) / 2;
+ }
+ }
+ _node_index->set(n, index);
+ (*_node_data)[index].pot = max;
+ int blossom =
+ _blossom_set->insert(n, std::numeric_limits::max());
+
+ _tree_set->insert(blossom);
+
+ (*_blossom_data)[blossom].status = EVEN;
+ (*_blossom_data)[blossom].pred = INVALID;
+ (*_blossom_data)[blossom].next = INVALID;
+ (*_blossom_data)[blossom].pot = 0;
+ (*_blossom_data)[blossom].offset = 0;
+ ++index;
+ }
+ for (EdgeIt e(_graph); e != INVALID; ++e) {
+ int si = (*_node_index)[_graph.u(e)];
+ int ti = (*_node_index)[_graph.v(e)];
+ if (_graph.u(e) != _graph.v(e)) {
+ _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
+ dualScale * _weight[e]) / 2);
+ }
+ }
+ }
+
+ /// \brief Starts the algorithm
+ ///
+ /// Starts the algorithm
+ bool start() {
+ enum OpType {
+ D2, D3, D4
+ };
+
+ int unmatched = _node_num;
+ while (unmatched > 0) {
+ Value d2 = !_delta2->empty() ?
+ _delta2->prio() : std::numeric_limits::max();
+
+ Value d3 = !_delta3->empty() ?
+ _delta3->prio() : std::numeric_limits::max();
+
+ Value d4 = !_delta4->empty() ?
+ _delta4->prio() : std::numeric_limits::max();
+
+ _delta_sum = d2; OpType ot = D2;
+ if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
+ if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
+
+ if (_delta_sum == std::numeric_limits::max()) {
+ return false;
+ }
+
+ switch (ot) {
+ case D2:
+ {
+ int blossom = _delta2->top();
+ Node n = _blossom_set->classTop(blossom);
+ Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
+ extendOnArc(e);
+ }
+ break;
+ case D3:
+ {
+ Edge e = _delta3->top();
+
+ int left_blossom = _blossom_set->find(_graph.u(e));
+ int right_blossom = _blossom_set->find(_graph.v(e));
+
+ if (left_blossom == right_blossom) {
+ _delta3->pop();
+ } else {
+ int left_tree = _tree_set->find(left_blossom);
+ int right_tree = _tree_set->find(right_blossom);
+
+ if (left_tree == right_tree) {
+ shrinkOnEdge(e, left_tree);
+ } else {
+ augmentOnEdge(e);
+ unmatched -= 2;
+ }
+ }
+ } break;
+ case D4:
+ splitBlossom(_delta4->top());
+ break;
+ }
+ }
+ extractMatching();
+ return true;
+ }
+
+ /// \brief Runs %MaxWeightedPerfectMatching algorithm.
+ ///
+ /// This method runs the %MaxWeightedPerfectMatching algorithm.
+ ///
+ /// \note mwm.run() is just a shortcut of the following code.
+ /// \code
+ /// mwm.init();
+ /// mwm.start();
+ /// \endcode
+ bool run() {
+ init();
+ return start();
+ }
+
+ /// @}
+
+ /// \name Primal solution
+ /// Functions to get the primal solution, ie. the matching.
+
+ /// @{
+
+ /// \brief Returns the matching value.
+ ///
+ /// Returns the matching value.
+ Value matchingValue() const {
+ Value sum = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if ((*_matching)[n] != INVALID) {
+ sum += _weight[(*_matching)[n]];
+ }
+ }
+ return sum /= 2;
+ }
+
+ /// \brief Returns true when the edge is in the matching.
+ ///
+ /// Returns true when the edge is in the matching.
+ bool matching(const Edge& edge) const {
+ return static_cast((*_matching)[_graph.u(edge)]) == edge;
+ }
+
+ /// \brief Returns the incident matching edge.
+ ///
+ /// Returns the incident matching arc from given edge.
+ Arc matching(const Node& node) const {
+ return (*_matching)[node];
+ }
+
+ /// \brief Returns the mate of the node.
+ ///
+ /// Returns the adjancent node in a mathcing arc.
+ Node mate(const Node& node) const {
+ return _graph.target((*_matching)[node]);
+ }
+
+ /// @}
+
+ /// \name Dual solution
+ /// Functions to get the dual solution.
+
+ /// @{
+
+ /// \brief Returns the value of the dual solution.
+ ///
+ /// Returns the value of the dual solution. It should be equal to
+ /// the primal value scaled by \ref dualScale "dual scale".
+ Value dualValue() const {
+ Value sum = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ sum += nodeValue(n);
+ }
+ for (int i = 0; i < blossomNum(); ++i) {
+ sum += blossomValue(i) * (blossomSize(i) / 2);
+ }
+ return sum;
+ }
+
+ /// \brief Returns the value of the node.
+ ///
+ /// Returns the the value of the node.
+ Value nodeValue(const Node& n) const {
+ return (*_node_potential)[n];
+ }
+
+ /// \brief Returns the number of the blossoms in the basis.
+ ///
+ /// Returns the number of the blossoms in the basis.
+ /// \see BlossomIt
+ int blossomNum() const {
+ return _blossom_potential.size();
+ }
+
+
+ /// \brief Returns the number of the nodes in the blossom.
+ ///
+ /// Returns the number of the nodes in the blossom.
+ int blossomSize(int k) const {
+ return _blossom_potential[k].end - _blossom_potential[k].begin;
+ }
+
+ /// \brief Returns the value of the blossom.
+ ///
+ /// Returns the the value of the blossom.
+ /// \see BlossomIt
+ Value blossomValue(int k) const {
+ return _blossom_potential[k].value;
+ }
+
+ /// \brief Iterator for obtaining the nodes of the blossom.
+ ///
+ /// Iterator for obtaining the nodes of the blossom. This class
+ /// provides a common lemon style iterator for listing a
+ /// subset of the nodes.
+ class BlossomIt {
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// Constructor to get the nodes of the variable.
+ BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
+ : _algorithm(&algorithm)
+ {
+ _index = _algorithm->_blossom_potential[variable].begin;
+ _last = _algorithm->_blossom_potential[variable].end;
+ }
+
+ /// \brief Conversion to node.
+ ///
+ /// Conversion to node.
+ operator Node() const {
+ return _algorithm->_blossom_node_list[_index];
+ }
+
+ /// \brief Increment operator.
+ ///
+ /// Increment operator.
+ BlossomIt& operator++() {
+ ++_index;
+ return *this;
+ }
+
+ /// \brief Validity checking
+ ///
+ /// Checks whether the iterator is invalid.
+ bool operator==(Invalid) const { return _index == _last; }
+
+ /// \brief Validity checking
+ ///
+ /// Checks whether the iterator is valid.
+ bool operator!=(Invalid) const { return _index != _last; }
+
+ private:
+ const MaxWeightedPerfectMatching* _algorithm;
+ int _last;
+ int _index;
+ };
+
+ /// @}
+
+ };
+
+
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_MAX_MATCHING_H
Index: lemon/nauty_reader.h
===================================================================
--- lemon/nauty_reader.h (revision 371)
+++ lemon/nauty_reader.h (revision 371)
@@ -0,0 +1,113 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_NAUTY_READER_H
+#define LEMON_NAUTY_READER_H
+
+#include
+#include
+#include
+
+/// \ingroup nauty_group
+/// \file
+/// \brief Nauty file reader.
+
+namespace lemon {
+
+ /// \ingroup nauty_group
+ ///
+ /// \brief Nauty file reader
+ ///
+ /// The \e geng program is in the \e gtools suite of the nauty
+ /// package. This tool can generate all non-isomorphic undirected
+ /// graphs of several classes with given node number (e.g.
+ /// general, connected, biconnected, triangle-free, 4-cycle-free,
+ /// bipartite and graphs with given edge number and degree
+ /// constraints). This function reads a \e nauty \e graph6 \e format
+ /// line from the given stream and builds it in the given graph.
+ ///
+ /// The site of nauty package: http://cs.anu.edu.au/~bdm/nauty/
+ ///
+ /// For example, the number of all non-isomorphic planar graphs
+ /// can be computed with the following code.
+ ///\code
+ /// int num = 0;
+ /// SmartGraph graph;
+ /// while (readNautyGraph(graph, std::cin)) {
+ /// PlanarityChecking pc(graph);
+ /// if (pc.run()) ++num;
+ /// }
+ /// std::cout << "Number of planar graphs: " << num << std::endl;
+ ///\endcode
+ ///
+ /// The nauty files are quite huge, therefore instead of the direct
+ /// file generation pipelining is recommended. For example,
+ ///\code
+ /// ./geng -c 10 | ./num_of_planar_graphs
+ ///\endcode
+ template
+ std::istream& readNautyGraph(Graph& graph, std::istream& is = std::cin) {
+ graph.clear();
+
+ std::string line;
+ if (getline(is, line)) {
+ int index = 0;
+
+ int n;
+
+ if (line[index] == '>') {
+ index += 10;
+ }
+
+ char c = line[index++]; c -= 63;
+ if (c != 63) {
+ n = int(c);
+ } else {
+ c = line[index++]; c -= 63;
+ n = (int(c) << 12);
+ c = line[index++]; c -= 63;
+ n |= (int(c) << 6);
+ c = line[index++]; c -= 63;
+ n |= int(c);
+ }
+
+ std::vector nodes;
+ for (int i = 0; i < n; ++i) {
+ nodes.push_back(graph.addNode());
+ }
+
+ int bit = -1;
+ for (int j = 0; j < n; ++j) {
+ for (int i = 0; i < j; ++i) {
+ if (bit == -1) {
+ c = line[index++]; c -= 63;
+ bit = 5;
+ }
+ bool b = (c & (1 << (bit--))) != 0;
+
+ if (b) {
+ graph.addEdge(nodes[i], nodes[j]);
+ }
+ }
+ }
+ }
+ return is;
+ }
+}
+
+#endif
Index: lemon/preflow.h
===================================================================
--- lemon/preflow.h (revision 408)
+++ lemon/preflow.h (revision 408)
@@ -0,0 +1,964 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_PREFLOW_H
+#define LEMON_PREFLOW_H
+
+#include
+#include
+
+/// \file
+/// \ingroup max_flow
+/// \brief Implementation of the preflow algorithm.
+
+namespace lemon {
+
+ /// \brief Default traits class of Preflow class.
+ ///
+ /// Default traits class of Preflow class.
+ /// \tparam _Digraph Digraph type.
+ /// \tparam _CapacityMap Capacity map type.
+ template
+ struct PreflowDefaultTraits {
+
+ /// \brief The type of the digraph the algorithm runs on.
+ typedef _Digraph Digraph;
+
+ /// \brief The type of the map that stores the arc capacities.
+ ///
+ /// The type of the map that stores the arc capacities.
+ /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+ typedef _CapacityMap CapacityMap;
+
+ /// \brief The type of the flow values.
+ typedef typename CapacityMap::Value Value;
+
+ /// \brief The type of the map that stores the flow values.
+ ///
+ /// The type of the map that stores the flow values.
+ /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+ typedef typename Digraph::template ArcMap FlowMap;
+
+ /// \brief Instantiates a FlowMap.
+ ///
+ /// This function instantiates a \ref FlowMap.
+ /// \param digraph The digraph, to which we would like to define
+ /// the flow map.
+ static FlowMap* createFlowMap(const Digraph& digraph) {
+ return new FlowMap(digraph);
+ }
+
+ /// \brief The elevator type used by Preflow algorithm.
+ ///
+ /// The elevator type used by Preflow algorithm.
+ ///
+ /// \sa Elevator
+ /// \sa LinkedElevator
+ typedef LinkedElevator Elevator;
+
+ /// \brief Instantiates an Elevator.
+ ///
+ /// This function instantiates an \ref Elevator.
+ /// \param digraph The digraph, to which we would like to define
+ /// the elevator.
+ /// \param max_level The maximum level of the elevator.
+ static Elevator* createElevator(const Digraph& digraph, int max_level) {
+ return new Elevator(digraph, max_level);
+ }
+
+ /// \brief The tolerance used by the algorithm
+ ///
+ /// The tolerance used by the algorithm to handle inexact computation.
+ typedef lemon::Tolerance Tolerance;
+
+ };
+
+
+ /// \ingroup max_flow
+ ///
+ /// \brief %Preflow algorithm class.
+ ///
+ /// This class provides an implementation of Goldberg-Tarjan's \e preflow
+ /// \e push-relabel algorithm producing a flow of maximum value in a
+ /// digraph. The preflow algorithms are the fastest known maximum
+ /// flow algorithms. The current implementation use a mixture of the
+ /// \e "highest label" and the \e "bound decrease" heuristics.
+ /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
+ ///
+ /// The algorithm consists of two phases. After the first phase
+ /// the maximum flow value and the minimum cut is obtained. The
+ /// second phase constructs a feasible maximum flow on each arc.
+ ///
+ /// \tparam _Digraph The type of the digraph the algorithm runs on.
+ /// \tparam _CapacityMap The type of the capacity map. The default map
+ /// type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap".
+#ifdef DOXYGEN
+ template
+#else
+ template ,
+ typename _Traits = PreflowDefaultTraits<_Digraph, _CapacityMap> >
+#endif
+ class Preflow {
+ public:
+
+ ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
+ typedef _Traits Traits;
+ ///The type of the digraph the algorithm runs on.
+ typedef typename Traits::Digraph Digraph;
+ ///The type of the capacity map.
+ typedef typename Traits::CapacityMap CapacityMap;
+ ///The type of the flow values.
+ typedef typename Traits::Value Value;
+
+ ///The type of the flow map.
+ typedef typename Traits::FlowMap FlowMap;
+ ///The type of the elevator.
+ typedef typename Traits::Elevator Elevator;
+ ///The type of the tolerance.
+ typedef typename Traits::Tolerance Tolerance;
+
+ private:
+
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+ const Digraph& _graph;
+ const CapacityMap* _capacity;
+
+ int _node_num;
+
+ Node _source, _target;
+
+ FlowMap* _flow;
+ bool _local_flow;
+
+ Elevator* _level;
+ bool _local_level;
+
+ typedef typename Digraph::template NodeMap ExcessMap;
+ ExcessMap* _excess;
+
+ Tolerance _tolerance;
+
+ bool _phase;
+
+
+ void createStructures() {
+ _node_num = countNodes(_graph);
+
+ if (!_flow) {
+ _flow = Traits::createFlowMap(_graph);
+ _local_flow = true;
+ }
+ if (!_level) {
+ _level = Traits::createElevator(_graph, _node_num);
+ _local_level = true;
+ }
+ if (!_excess) {
+ _excess = new ExcessMap(_graph);
+ }
+ }
+
+ void destroyStructures() {
+ if (_local_flow) {
+ delete _flow;
+ }
+ if (_local_level) {
+ delete _level;
+ }
+ if (_excess) {
+ delete _excess;
+ }
+ }
+
+ public:
+
+ typedef Preflow Create;
+
+ ///\name Named Template Parameters
+
+ ///@{
+
+ template
+ struct SetFlowMapTraits : public Traits {
+ typedef _FlowMap FlowMap;
+ static FlowMap *createFlowMap(const Digraph&) {
+ LEMON_ASSERT(false, "FlowMap is not initialized");
+ return 0; // ignore warnings
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// FlowMap type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting FlowMap
+ /// type.
+ template
+ struct SetFlowMap
+ : public Preflow > {
+ typedef Preflow > Create;
+ };
+
+ template
+ struct SetElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Digraph&, int) {
+ LEMON_ASSERT(false, "Elevator is not initialized");
+ return 0; // ignore warnings
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// Elevator type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting Elevator
+ /// type. If this named parameter is used, then an external
+ /// elevator object must be passed to the algorithm using the
+ /// \ref elevator(Elevator&) "elevator()" function before calling
+ /// \ref run() or \ref init().
+ /// \sa SetStandardElevator
+ template
+ struct SetElevator
+ : public Preflow > {
+ typedef Preflow > Create;
+ };
+
+ template
+ struct SetStandardElevatorTraits : public Traits {
+ typedef _Elevator Elevator;
+ static Elevator *createElevator(const Digraph& digraph, int max_level) {
+ return new Elevator(digraph, max_level);
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting
+ /// Elevator type with automatic allocation
+ ///
+ /// \ref named-templ-param "Named parameter" for setting Elevator
+ /// type with automatic allocation.
+ /// The Elevator should have standard constructor interface to be
+ /// able to automatically created by the algorithm (i.e. the
+ /// digraph and the maximum level should be passed to it).
+ /// However an external elevator object could also be passed to the
+ /// algorithm with the \ref elevator(Elevator&) "elevator()" function
+ /// before calling \ref run() or \ref init().
+ /// \sa SetElevator
+ template
+ struct SetStandardElevator
+ : public Preflow > {
+ typedef Preflow > Create;
+ };
+
+ /// @}
+
+ protected:
+
+ Preflow() {}
+
+ public:
+
+
+ /// \brief The constructor of the class.
+ ///
+ /// The constructor of the class.
+ /// \param digraph The digraph the algorithm runs on.
+ /// \param capacity The capacity of the arcs.
+ /// \param source The source node.
+ /// \param target The target node.
+ Preflow(const Digraph& digraph, const CapacityMap& capacity,
+ Node source, Node target)
+ : _graph(digraph), _capacity(&capacity),
+ _node_num(0), _source(source), _target(target),
+ _flow(0), _local_flow(false),
+ _level(0), _local_level(false),
+ _excess(0), _tolerance(), _phase() {}
+
+ /// \brief Destructor.
+ ///
+ /// Destructor.
+ ~Preflow() {
+ destroyStructures();
+ }
+
+ /// \brief Sets the capacity map.
+ ///
+ /// Sets the capacity map.
+ /// \return (*this)
+ Preflow& capacityMap(const CapacityMap& map) {
+ _capacity = ↦
+ return *this;
+ }
+
+ /// \brief Sets the flow map.
+ ///
+ /// Sets the flow map.
+ /// If you don't use this function before calling \ref run() or
+ /// \ref init(), an instance will be allocated automatically.
+ /// The destructor deallocates this automatically allocated map,
+ /// of course.
+ /// \return (*this)
+ Preflow& flowMap(FlowMap& map) {
+ if (_local_flow) {
+ delete _flow;
+ _local_flow = false;
+ }
+ _flow = ↦
+ return *this;
+ }
+
+ /// \brief Sets the source node.
+ ///
+ /// Sets the source node.
+ /// \return (*this)
+ Preflow& source(const Node& node) {
+ _source = node;
+ return *this;
+ }
+
+ /// \brief Sets the target node.
+ ///
+ /// Sets the target node.
+ /// \return (*this)
+ Preflow& target(const Node& node) {
+ _target = node;
+ return *this;
+ }
+
+ /// \brief Sets the elevator used by algorithm.
+ ///
+ /// Sets the elevator used by algorithm.
+ /// If you don't use this function before calling \ref run() or
+ /// \ref init(), an instance will be allocated automatically.
+ /// The destructor deallocates this automatically allocated elevator,
+ /// of course.
+ /// \return (*this)
+ Preflow& elevator(Elevator& elevator) {
+ if (_local_level) {
+ delete _level;
+ _local_level = false;
+ }
+ _level = &elevator;
+ return *this;
+ }
+
+ /// \brief Returns a const reference to the elevator.
+ ///
+ /// Returns a const reference to the elevator.
+ ///
+ /// \pre Either \ref run() or \ref init() must be called before
+ /// using this function.
+ const Elevator& elevator() {
+ return *_level;
+ }
+
+ /// \brief Sets the tolerance used by algorithm.
+ ///
+ /// Sets the tolerance used by algorithm.
+ Preflow& tolerance(const Tolerance& tolerance) const {
+ _tolerance = tolerance;
+ return *this;
+ }
+
+ /// \brief Returns a const reference to the tolerance.
+ ///
+ /// Returns a const reference to the tolerance.
+ const Tolerance& tolerance() const {
+ return tolerance;
+ }
+
+ /// \name Execution Control
+ /// The simplest way to execute the preflow algorithm is to use
+ /// \ref run() or \ref runMinCut().\n
+ /// If you need more control on the initial solution or the execution,
+ /// first you have to call one of the \ref init() functions, then
+ /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
+
+ ///@{
+
+ /// \brief Initializes the internal data structures.
+ ///
+ /// Initializes the internal data structures and sets the initial
+ /// flow to zero on each arc.
+ void init() {
+ createStructures();
+
+ _phase = true;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _excess->set(n, 0);
+ }
+
+ for (ArcIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, 0);
+ }
+
+ typename Digraph::template NodeMap reached(_graph, false);
+
+ _level->initStart();
+ _level->initAddItem(_target);
+
+ std::vector queue;
+ reached.set(_source, true);
+
+ queue.push_back(_target);
+ reached.set(_target, true);
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
+ }
+ }
+ }
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
+
+ for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
+ if (_tolerance.positive((*_capacity)[e])) {
+ Node u = _graph.target(e);
+ if ((*_level)[u] == _level->maxLevel()) continue;
+ _flow->set(e, (*_capacity)[e]);
+ _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+ if (u != _target && !_level->active(u)) {
+ _level->activate(u);
+ }
+ }
+ }
+ }
+
+ /// \brief Initializes the internal data structures using the
+ /// given flow map.
+ ///
+ /// Initializes the internal data structures and sets the initial
+ /// flow to the given \c flowMap. The \c flowMap should contain a
+ /// flow or at least a preflow, i.e. at each node excluding the
+ /// source node the incoming flow should greater or equal to the
+ /// outgoing flow.
+ /// \return \c false if the given \c flowMap is not a preflow.
+ template
+ bool init(const FlowMap& flowMap) {
+ createStructures();
+
+ for (ArcIt e(_graph); e != INVALID; ++e) {
+ _flow->set(e, flowMap[e]);
+ }
+
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ Value excess = 0;
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ excess += (*_flow)[e];
+ }
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ excess -= (*_flow)[e];
+ }
+ if (excess < 0 && n != _source) return false;
+ _excess->set(n, excess);
+ }
+
+ typename Digraph::template NodeMap reached(_graph, false);
+
+ _level->initStart();
+ _level->initAddItem(_target);
+
+ std::vector queue;
+ reached.set(_source, true);
+
+ queue.push_back(_target);
+ reached.set(_target, true);
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] &&
+ _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
+ }
+ }
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+ reached.set(v, true);
+ _level->initAddItem(v);
+ nqueue.push_back(v);
+ }
+ }
+ }
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
+
+ for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (_tolerance.positive(rem)) {
+ Node u = _graph.target(e);
+ if ((*_level)[u] == _level->maxLevel()) continue;
+ _flow->set(e, (*_capacity)[e]);
+ _excess->set(u, (*_excess)[u] + rem);
+ if (u != _target && !_level->active(u)) {
+ _level->activate(u);
+ }
+ }
+ }
+ for (InArcIt e(_graph, _source); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (_tolerance.positive(rem)) {
+ Node v = _graph.source(e);
+ if ((*_level)[v] == _level->maxLevel()) continue;
+ _flow->set(e, 0);
+ _excess->set(v, (*_excess)[v] + rem);
+ if (v != _target && !_level->active(v)) {
+ _level->activate(v);
+ }
+ }
+ }
+ return true;
+ }
+
+ /// \brief Starts the first phase of the preflow algorithm.
+ ///
+ /// The preflow algorithm consists of two phases, this method runs
+ /// the first phase. After the first phase the maximum flow value
+ /// and a minimum value cut can already be computed, although a
+ /// maximum flow is not yet obtained. So after calling this method
+ /// \ref flowValue() returns the value of a maximum flow and \ref
+ /// minCut() returns a minimum cut.
+ /// \pre One of the \ref init() functions must be called before
+ /// using this function.
+ void startFirstPhase() {
+ _phase = true;
+
+ Node n = _level->highestActive();
+ int level = _level->highestActiveLevel();
+ while (n != INVALID) {
+ int num = _node_num;
+
+ while (num > 0 && n != INVALID) {
+ Value excess = (*_excess)[n];
+ int new_level = _level->maxLevel();
+
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.target(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_1;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.source(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_1;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ no_more_push_1:
+
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftHighestActive(new_level + 1);
+ } else {
+ _level->liftHighestActiveToTop();
+ }
+ if (_level->emptyLevel(level)) {
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
+
+ n = _level->highestActive();
+ level = _level->highestActiveLevel();
+ --num;
+ }
+
+ num = _node_num * 20;
+ while (num > 0 && n != INVALID) {
+ Value excess = (*_excess)[n];
+ int new_level = _level->maxLevel();
+
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.target(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_2;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.source(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _target) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push_2;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ no_more_push_2:
+
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftActiveOn(level, new_level + 1);
+ } else {
+ _level->liftActiveToTop(level);
+ }
+ if (_level->emptyLevel(level)) {
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
+
+ while (level >= 0 && _level->activeFree(level)) {
+ --level;
+ }
+ if (level == -1) {
+ n = _level->highestActive();
+ level = _level->highestActiveLevel();
+ } else {
+ n = _level->activeOn(level);
+ }
+ --num;
+ }
+ }
+ }
+
+ /// \brief Starts the second phase of the preflow algorithm.
+ ///
+ /// The preflow algorithm consists of two phases, this method runs
+ /// the second phase. After calling one of the \ref init() functions
+ /// and \ref startFirstPhase() and then \ref startSecondPhase(),
+ /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
+ /// value of a maximum flow, \ref minCut() returns a minimum cut
+ /// \pre One of the \ref init() functions and \ref startFirstPhase()
+ /// must be called before using this function.
+ void startSecondPhase() {
+ _phase = false;
+
+ typename Digraph::template NodeMap reached(_graph);
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ reached.set(n, (*_level)[n] < _level->maxLevel());
+ }
+
+ _level->initStart();
+ _level->initAddItem(_source);
+
+ std::vector queue;
+ queue.push_back(_source);
+ reached.set(_source, true);
+
+ while (!queue.empty()) {
+ _level->initNewLevel();
+ std::vector nqueue;
+ for (int i = 0; i < int(queue.size()); ++i) {
+ Node n = queue[i];
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ Node v = _graph.target(e);
+ if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+ reached.set(v, true);
+ _level->initAddItem(v);
+ nqueue.push_back(v);
+ }
+ }
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Node u = _graph.source(e);
+ if (!reached[u] &&
+ _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+ reached.set(u, true);
+ _level->initAddItem(u);
+ nqueue.push_back(u);
+ }
+ }
+ }
+ queue.swap(nqueue);
+ }
+ _level->initFinish();
+
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if (!reached[n]) {
+ _level->dirtyTopButOne(n);
+ } else if ((*_excess)[n] > 0 && _target != n) {
+ _level->activate(n);
+ }
+ }
+
+ Node n;
+ while ((n = _level->highestActive()) != INVALID) {
+ Value excess = (*_excess)[n];
+ int level = _level->highestActiveLevel();
+ int new_level = _level->maxLevel();
+
+ for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_capacity)[e] - (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.target(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _source) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] + excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, (*_capacity)[e]);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ for (InArcIt e(_graph, n); e != INVALID; ++e) {
+ Value rem = (*_flow)[e];
+ if (!_tolerance.positive(rem)) continue;
+ Node v = _graph.source(e);
+ if ((*_level)[v] < level) {
+ if (!_level->active(v) && v != _source) {
+ _level->activate(v);
+ }
+ if (!_tolerance.less(rem, excess)) {
+ _flow->set(e, (*_flow)[e] - excess);
+ _excess->set(v, (*_excess)[v] + excess);
+ excess = 0;
+ goto no_more_push;
+ } else {
+ excess -= rem;
+ _excess->set(v, (*_excess)[v] + rem);
+ _flow->set(e, 0);
+ }
+ } else if (new_level > (*_level)[v]) {
+ new_level = (*_level)[v];
+ }
+ }
+
+ no_more_push:
+
+ _excess->set(n, excess);
+
+ if (excess != 0) {
+ if (new_level + 1 < _level->maxLevel()) {
+ _level->liftHighestActive(new_level + 1);
+ } else {
+ // Calculation error
+ _level->liftHighestActiveToTop();
+ }
+ if (_level->emptyLevel(level)) {
+ // Calculation error
+ _level->liftToTop(level);
+ }
+ } else {
+ _level->deactivate(n);
+ }
+
+ }
+ }
+
+ /// \brief Runs the preflow algorithm.
+ ///
+ /// Runs the preflow algorithm.
+ /// \note pf.run() is just a shortcut of the following code.
+ /// \code
+ /// pf.init();
+ /// pf.startFirstPhase();
+ /// pf.startSecondPhase();
+ /// \endcode
+ void run() {
+ init();
+ startFirstPhase();
+ startSecondPhase();
+ }
+
+ /// \brief Runs the preflow algorithm to compute the minimum cut.
+ ///
+ /// Runs the preflow algorithm to compute the minimum cut.
+ /// \note pf.runMinCut() is just a shortcut of the following code.
+ /// \code
+ /// pf.init();
+ /// pf.startFirstPhase();
+ /// \endcode
+ void runMinCut() {
+ init();
+ startFirstPhase();
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The results of the preflow algorithm can be obtained using these
+ /// functions.\n
+ /// Either one of the \ref run() "run*()" functions or one of the
+ /// \ref startFirstPhase() "start*()" functions should be called
+ /// before using them.
+
+ ///@{
+
+ /// \brief Returns the value of the maximum flow.
+ ///
+ /// Returns the value of the maximum flow by returning the excess
+ /// of the target node. This value equals to the value of
+ /// the maximum flow already after the first phase of the algorithm.
+ ///
+ /// \pre Either \ref run() or \ref init() must be called before
+ /// using this function.
+ Value flowValue() const {
+ return (*_excess)[_target];
+ }
+
+ /// \brief Returns the flow on the given arc.
+ ///
+ /// Returns the flow on the given arc. This method can
+ /// be called after the second phase of the algorithm.
+ ///
+ /// \pre Either \ref run() or \ref init() must be called before
+ /// using this function.
+ Value flow(const Arc& arc) const {
+ return (*_flow)[arc];
+ }
+
+ /// \brief Returns a const reference to the flow map.
+ ///
+ /// Returns a const reference to the arc map storing the found flow.
+ /// This method can be called after the second phase of the algorithm.
+ ///
+ /// \pre Either \ref run() or \ref init() must be called before
+ /// using this function.
+ const FlowMap& flowMap() {
+ return *_flow;
+ }
+
+ /// \brief Returns \c true when the node is on the source side of the
+ /// minimum cut.
+ ///
+ /// Returns true when the node is on the source side of the found
+ /// minimum cut. This method can be called both after running \ref
+ /// startFirstPhase() and \ref startSecondPhase().
+ ///
+ /// \pre Either \ref run() or \ref init() must be called before
+ /// using this function.
+ bool minCut(const Node& node) const {
+ return ((*_level)[node] == _level->maxLevel()) == _phase;
+ }
+
+ /// \brief Gives back a minimum value cut.
+ ///
+ /// Sets \c cutMap to the characteristic vector of a minimum value
+ /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
+ /// node map with \c bool (or convertible) value type.
+ ///
+ /// This method can be called both after running \ref startFirstPhase()
+ /// and \ref startSecondPhase(). The result after the second phase
+ /// could be slightly different if inexact computation is used.
+ ///
+ /// \note This function calls \ref minCut() for each node, so it runs in
+ /// \f$O(n)\f$ time.
+ ///
+ /// \pre Either \ref run() or \ref init() must be called before
+ /// using this function.
+ template
+ void minCutMap(CutMap& cutMap) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ cutMap.set(n, minCut(n));
+ }
+ }
+
+ /// @}
+ };
+}
+
+#endif
Index: lemon/random.h
===================================================================
--- lemon/random.h (revision 391)
+++ lemon/random.h (revision 392)
@@ -541,8 +541,4 @@
/// @{
- ///\name Initialization
- ///
- /// @{
-
/// \brief Default constructor
///
@@ -693,10 +689,4 @@
}
- /// @}
-
- ///\name Uniform distributions
- ///
- /// @{
-
/// \brief Returns a random real number from the range [0, 1)
///
@@ -753,6 +743,4 @@
return _random_bits::IntConversion::convert(core);
}
-
- /// @}
unsigned int uinteger() {
@@ -789,8 +777,7 @@
///\name Non-uniform distributions
///
-
///@{
- /// \brief Returns a random bool
+ /// \brief Returns a random bool with given probability of true result.
///
/// It returns a random bool with given probability of true result.
@@ -799,7 +786,7 @@
}
- /// Standard Gauss distribution
-
- /// Standard Gauss distribution.
+ /// Standard normal (Gauss) distribution
+
+ /// Standard normal (Gauss) distribution.
/// \note The Cartesian form of the Box-Muller
/// transformation is used to generate a random normal distribution.
@@ -814,11 +801,51 @@
return std::sqrt(-2*std::log(S)/S)*V1;
}
- /// Gauss distribution with given mean and standard deviation
-
- /// Gauss distribution with given mean and standard deviation.
+ /// Normal (Gauss) distribution with given mean and standard deviation
+
+ /// Normal (Gauss) distribution with given mean and standard deviation.
/// \sa gauss()
double gauss(double mean,double std_dev)
{
return gauss()*std_dev+mean;
+ }
+
+ /// Lognormal distribution
+
+ /// Lognormal distribution. The parameters are the mean and the standard
+ /// deviation of exp(X).
+ ///
+ double lognormal(double n_mean,double n_std_dev)
+ {
+ return std::exp(gauss(n_mean,n_std_dev));
+ }
+ /// Lognormal distribution
+
+ /// Lognormal distribution. The parameter is an std::pair of
+ /// the mean and the standard deviation of exp(X).
+ ///
+ double lognormal(const std::pair ¶ms)
+ {
+ return std::exp(gauss(params.first,params.second));
+ }
+ /// Compute the lognormal parameters from mean and standard deviation
+
+ /// This function computes the lognormal parameters from mean and
+ /// standard deviation. The return value can direcly be passed to
+ /// lognormal().
+ std::pair lognormalParamsFromMD(double mean,
+ double std_dev)
+ {
+ double fr=std_dev/mean;
+ fr*=fr;
+ double lg=std::log(1+fr);
+ return std::pair(std::log(mean)-lg/2.0,std::sqrt(lg));
+ }
+ /// Lognormal distribution with given mean and standard deviation
+
+ /// Lognormal distribution with given mean and standard deviation.
+ ///
+ double lognormalMD(double mean,double std_dev)
+ {
+ return lognormal(lognormalParamsFromMD(mean,std_dev));
}
@@ -926,5 +953,4 @@
///\name Two dimensional distributions
///
-
///@{
@@ -943,5 +969,5 @@
return dim2::Point(V1,V2);
}
- /// A kind of two dimensional Gauss distribution
+ /// A kind of two dimensional normal (Gauss) distribution
/// This function provides a turning symmetric two-dimensional distribution.
Index: lemon/smart_graph.h
===================================================================
--- lemon/smart_graph.h (revision 313)
+++ lemon/smart_graph.h (revision 388)
@@ -68,5 +68,5 @@
typedef True NodeNumTag;
- typedef True EdgeNumTag;
+ typedef True ArcNumTag;
int nodeNum() const { return nodes.size(); }
@@ -306,5 +306,7 @@
nodes[b._id].first_out=nodes[n._id].first_out;
nodes[n._id].first_out=-1;
- for(int i=nodes[b._id].first_out;i!=-1;i++) arcs[i].source=b._id;
+ for(int i=nodes[b._id].first_out; i!=-1; i=arcs[i].next_out) {
+ arcs[i].source=b._id;
+ }
if(connect) addArc(n,b);
return b;
@@ -465,6 +467,6 @@
public:
- operator Edge() const {
- return _id != -1 ? edgeFromId(_id / 2) : INVALID;
+ operator Edge() const {
+ return _id != -1 ? edgeFromId(_id / 2) : INVALID;
}
@@ -481,4 +483,11 @@
: nodes(), arcs() {}
+ typedef True NodeNumTag;
+ typedef True EdgeNumTag;
+ typedef True ArcNumTag;
+
+ int nodeNum() const { return nodes.size(); }
+ int edgeNum() const { return arcs.size() / 2; }
+ int arcNum() const { return arcs.size(); }
int maxNodeId() const { return nodes.size()-1; }
@@ -729,6 +738,6 @@
dir.push_back(arcFromId(n-1));
Parent::notifier(Arc()).erase(dir);
- nodes[arcs[n].target].first_out=arcs[n].next_out;
- nodes[arcs[n-1].target].first_out=arcs[n-1].next_out;
+ nodes[arcs[n-1].target].first_out=arcs[n].next_out;
+ nodes[arcs[n].target].first_out=arcs[n-1].next_out;
arcs.pop_back();
arcs.pop_back();
Index: lemon/suurballe.h
===================================================================
--- lemon/suurballe.h (revision 358)
+++ lemon/suurballe.h (revision 358)
@@ -0,0 +1,501 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_SUURBALLE_H
+#define LEMON_SUURBALLE_H
+
+///\ingroup shortest_path
+///\file
+///\brief An algorithm for finding arc-disjoint paths between two
+/// nodes having minimum total length.
+
+#include
+#include
+#include
+
+namespace lemon {
+
+ /// \addtogroup shortest_path
+ /// @{
+
+ /// \brief Algorithm for finding arc-disjoint paths between two nodes
+ /// having minimum total length.
+ ///
+ /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
+ /// finding arc-disjoint paths having minimum total length (cost)
+ /// from a given source node to a given target node in a digraph.
+ ///
+ /// In fact, this implementation is the specialization of the
+ /// \ref CapacityScaling "successive shortest path" algorithm.
+ ///
+ /// \tparam Digraph The digraph type the algorithm runs on.
+ /// The default value is \c ListDigraph.
+ /// \tparam LengthMap The type of the length (cost) map.
+ /// The default value is Digraph::ArcMap.
+ ///
+ /// \warning Length values should be \e non-negative \e integers.
+ ///
+ /// \note For finding node-disjoint paths this algorithm can be used
+ /// with \ref SplitDigraphAdaptor.
+#ifdef DOXYGEN
+ template
+#else
+ template < typename Digraph = ListDigraph,
+ typename LengthMap = typename Digraph::template ArcMap >
+#endif
+ class Suurballe
+ {
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+ typedef typename LengthMap::Value Length;
+ typedef ConstMap ConstArcMap;
+ typedef typename Digraph::template NodeMap PredMap;
+
+ public:
+
+ /// The type of the flow map.
+ typedef typename Digraph::template ArcMap FlowMap;
+ /// The type of the potential map.
+ typedef typename Digraph::template NodeMap PotentialMap;
+ /// The type of the path structures.
+ typedef SimplePath Path;
+
+ private:
+
+ /// \brief Special implementation of the Dijkstra algorithm
+ /// for finding shortest paths in the residual network.
+ ///
+ /// \ref ResidualDijkstra is a special implementation of the
+ /// \ref Dijkstra algorithm for finding shortest paths in the
+ /// residual network of the digraph with respect to the reduced arc
+ /// lengths and modifying the node potentials according to the
+ /// distance of the nodes.
+ class ResidualDijkstra
+ {
+ typedef typename Digraph::template NodeMap HeapCrossRef;
+ typedef BinHeap Heap;
+
+ private:
+
+ // The digraph the algorithm runs on
+ const Digraph &_graph;
+
+ // The main maps
+ const FlowMap &_flow;
+ const LengthMap &_length;
+ PotentialMap &_potential;
+
+ // The distance map
+ PotentialMap _dist;
+ // The pred arc map
+ PredMap &_pred;
+ // The processed (i.e. permanently labeled) nodes
+ std::vector _proc_nodes;
+
+ Node _s;
+ Node _t;
+
+ public:
+
+ /// Constructor.
+ ResidualDijkstra( const Digraph &digraph,
+ const FlowMap &flow,
+ const LengthMap &length,
+ PotentialMap &potential,
+ PredMap &pred,
+ Node s, Node t ) :
+ _graph(digraph), _flow(flow), _length(length), _potential(potential),
+ _dist(digraph), _pred(pred), _s(s), _t(t) {}
+
+ /// \brief Run the algorithm. It returns \c true if a path is found
+ /// from the source node to the target node.
+ bool run() {
+ HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
+ Heap heap(heap_cross_ref);
+ heap.push(_s, 0);
+ _pred[_s] = INVALID;
+ _proc_nodes.clear();
+
+ // Process nodes
+ while (!heap.empty() && heap.top() != _t) {
+ Node u = heap.top(), v;
+ Length d = heap.prio() + _potential[u], nd;
+ _dist[u] = heap.prio();
+ heap.pop();
+ _proc_nodes.push_back(u);
+
+ // Traverse outgoing arcs
+ for (OutArcIt e(_graph, u); e != INVALID; ++e) {
+ if (_flow[e] == 0) {
+ v = _graph.target(e);
+ switch(heap.state(v)) {
+ case Heap::PRE_HEAP:
+ heap.push(v, d + _length[e] - _potential[v]);
+ _pred[v] = e;
+ break;
+ case Heap::IN_HEAP:
+ nd = d + _length[e] - _potential[v];
+ if (nd < heap[v]) {
+ heap.decrease(v, nd);
+ _pred[v] = e;
+ }
+ break;
+ case Heap::POST_HEAP:
+ break;
+ }
+ }
+ }
+
+ // Traverse incoming arcs
+ for (InArcIt e(_graph, u); e != INVALID; ++e) {
+ if (_flow[e] == 1) {
+ v = _graph.source(e);
+ switch(heap.state(v)) {
+ case Heap::PRE_HEAP:
+ heap.push(v, d - _length[e] - _potential[v]);
+ _pred[v] = e;
+ break;
+ case Heap::IN_HEAP:
+ nd = d - _length[e] - _potential[v];
+ if (nd < heap[v]) {
+ heap.decrease(v, nd);
+ _pred[v] = e;
+ }
+ break;
+ case Heap::POST_HEAP:
+ break;
+ }
+ }
+ }
+ }
+ if (heap.empty()) return false;
+
+ // Update potentials of processed nodes
+ Length t_dist = heap.prio();
+ for (int i = 0; i < int(_proc_nodes.size()); ++i)
+ _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
+ return true;
+ }
+
+ }; //class ResidualDijkstra
+
+ private:
+
+ // The digraph the algorithm runs on
+ const Digraph &_graph;
+ // The length map
+ const LengthMap &_length;
+
+ // Arc map of the current flow
+ FlowMap *_flow;
+ bool _local_flow;
+ // Node map of the current potentials
+ PotentialMap *_potential;
+ bool _local_potential;
+
+ // The source node
+ Node _source;
+ // The target node
+ Node _target;
+
+ // Container to store the found paths
+ std::vector< SimplePath > paths;
+ int _path_num;
+
+ // The pred arc map
+ PredMap _pred;
+ // Implementation of the Dijkstra algorithm for finding augmenting
+ // shortest paths in the residual network
+ ResidualDijkstra *_dijkstra;
+
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// Constructor.
+ ///
+ /// \param digraph The digraph the algorithm runs on.
+ /// \param length The length (cost) values of the arcs.
+ /// \param s The source node.
+ /// \param t The target node.
+ Suurballe( const Digraph &digraph,
+ const LengthMap &length,
+ Node s, Node t ) :
+ _graph(digraph), _length(length), _flow(0), _local_flow(false),
+ _potential(0), _local_potential(false), _source(s), _target(t),
+ _pred(digraph) {}
+
+ /// Destructor.
+ ~Suurballe() {
+ if (_local_flow) delete _flow;
+ if (_local_potential) delete _potential;
+ delete _dijkstra;
+ }
+
+ /// \brief Set the flow map.
+ ///
+ /// This function sets the flow map.
+ ///
+ /// The found flow contains only 0 and 1 values. It is the union of
+ /// the found arc-disjoint paths.
+ ///
+ /// \return \c (*this)
+ Suurballe& flowMap(FlowMap &map) {
+ if (_local_flow) {
+ delete _flow;
+ _local_flow = false;
+ }
+ _flow = ↦
+ return *this;
+ }
+
+ /// \brief Set the potential map.
+ ///
+ /// This function sets the potential map.
+ ///
+ /// The potentials provide the dual solution of the underlying
+ /// minimum cost flow problem.
+ ///
+ /// \return \c (*this)
+ Suurballe& potentialMap(PotentialMap &map) {
+ if (_local_potential) {
+ delete _potential;
+ _local_potential = false;
+ }
+ _potential = ↦
+ return *this;
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to call the run()
+ /// function.
+ /// \n
+ /// If you only need the flow that is the union of the found
+ /// arc-disjoint paths, you may call init() and findFlow().
+
+ /// @{
+
+ /// \brief Run the algorithm.
+ ///
+ /// This function runs the algorithm.
+ ///
+ /// \param k The number of paths to be found.
+ ///
+ /// \return \c k if there are at least \c k arc-disjoint paths from
+ /// \c s to \c t in the digraph. Otherwise it returns the number of
+ /// arc-disjoint paths found.
+ ///
+ /// \note Apart from the return value, s.run(k) is just a
+ /// shortcut of the following code.
+ /// \code
+ /// s.init();
+ /// s.findFlow(k);
+ /// s.findPaths();
+ /// \endcode
+ int run(int k = 2) {
+ init();
+ findFlow(k);
+ findPaths();
+ return _path_num;
+ }
+
+ /// \brief Initialize the algorithm.
+ ///
+ /// This function initializes the algorithm.
+ void init() {
+ // Initialize maps
+ if (!_flow) {
+ _flow = new FlowMap(_graph);
+ _local_flow = true;
+ }
+ if (!_potential) {
+ _potential = new PotentialMap(_graph);
+ _local_potential = true;
+ }
+ for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
+
+ _dijkstra = new ResidualDijkstra( _graph, *_flow, _length,
+ *_potential, _pred,
+ _source, _target );
+ }
+
+ /// \brief Execute the successive shortest path algorithm to find
+ /// an optimal flow.
+ ///
+ /// This function executes the successive shortest path algorithm to
+ /// find a minimum cost flow, which is the union of \c k or less
+ /// arc-disjoint paths.
+ ///
+ /// \return \c k if there are at least \c k arc-disjoint paths from
+ /// \c s to \c t in the digraph. Otherwise it returns the number of
+ /// arc-disjoint paths found.
+ ///
+ /// \pre \ref init() must be called before using this function.
+ int findFlow(int k = 2) {
+ // Find shortest paths
+ _path_num = 0;
+ while (_path_num < k) {
+ // Run Dijkstra
+ if (!_dijkstra->run()) break;
+ ++_path_num;
+
+ // Set the flow along the found shortest path
+ Node u = _target;
+ Arc e;
+ while ((e = _pred[u]) != INVALID) {
+ if (u == _graph.target(e)) {
+ (*_flow)[e] = 1;
+ u = _graph.source(e);
+ } else {
+ (*_flow)[e] = 0;
+ u = _graph.target(e);
+ }
+ }
+ }
+ return _path_num;
+ }
+
+ /// \brief Compute the paths from the flow.
+ ///
+ /// This function computes the paths from the flow.
+ ///
+ /// \pre \ref init() and \ref findFlow() must be called before using
+ /// this function.
+ void findPaths() {
+ // Create the residual flow map (the union of the paths not found
+ // so far)
+ FlowMap res_flow(_graph);
+ for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
+
+ paths.clear();
+ paths.resize(_path_num);
+ for (int i = 0; i < _path_num; ++i) {
+ Node n = _source;
+ while (n != _target) {
+ OutArcIt e(_graph, n);
+ for ( ; res_flow[e] == 0; ++e) ;
+ n = _graph.target(e);
+ paths[i].addBack(e);
+ res_flow[e] = 0;
+ }
+ }
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The results of the algorithm can be obtained using these
+ /// functions.
+ /// \n The algorithm should be executed before using them.
+
+ /// @{
+
+ /// \brief Return a const reference to the arc map storing the
+ /// found flow.
+ ///
+ /// This function returns a const reference to the arc map storing
+ /// the flow that is the union of the found arc-disjoint paths.
+ ///
+ /// \pre \ref run() or \ref findFlow() must be called before using
+ /// this function.
+ const FlowMap& flowMap() const {
+ return *_flow;
+ }
+
+ /// \brief Return a const reference to the node map storing the
+ /// found potentials (the dual solution).
+ ///
+ /// This function returns a const reference to the node map storing
+ /// the found potentials that provide the dual solution of the
+ /// underlying minimum cost flow problem.
+ ///
+ /// \pre \ref run() or \ref findFlow() must be called before using
+ /// this function.
+ const PotentialMap& potentialMap() const {
+ return *_potential;
+ }
+
+ /// \brief Return the flow on the given arc.
+ ///
+ /// This function returns the flow on the given arc.
+ /// It is \c 1 if the arc is involved in one of the found paths,
+ /// otherwise it is \c 0.
+ ///
+ /// \pre \ref run() or \ref findFlow() must be called before using
+ /// this function.
+ int flow(const Arc& arc) const {
+ return (*_flow)[arc];
+ }
+
+ /// \brief Return the potential of the given node.
+ ///
+ /// This function returns the potential of the given node.
+ ///
+ /// \pre \ref run() or \ref findFlow() must be called before using
+ /// this function.
+ Length potential(const Node& node) const {
+ return (*_potential)[node];
+ }
+
+ /// \brief Return the total length (cost) of the found paths (flow).
+ ///
+ /// This function returns the total length (cost) of the found paths
+ /// (flow). The complexity of the function is \f$ O(e) \f$.
+ ///
+ /// \pre \ref run() or \ref findFlow() must be called before using
+ /// this function.
+ Length totalLength() const {
+ Length c = 0;
+ for (ArcIt e(_graph); e != INVALID; ++e)
+ c += (*_flow)[e] * _length[e];
+ return c;
+ }
+
+ /// \brief Return the number of the found paths.
+ ///
+ /// This function returns the number of the found paths.
+ ///
+ /// \pre \ref run() or \ref findFlow() must be called before using
+ /// this function.
+ int pathNum() const {
+ return _path_num;
+ }
+
+ /// \brief Return a const reference to the specified path.
+ ///
+ /// This function returns a const reference to the specified path.
+ ///
+ /// \param i The function returns the \c i-th path.
+ /// \c i must be between \c 0 and %pathNum()-1.
+ ///
+ /// \pre \ref run() or \ref findPaths() must be called before using
+ /// this function.
+ Path path(int i) const {
+ return paths[i];
+ }
+
+ /// @}
+
+ }; //class Suurballe
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_SUURBALLE_H
Index: scripts/chg-len.py
===================================================================
--- scripts/chg-len.py (revision 284)
+++ scripts/chg-len.py (revision 390)
@@ -2,5 +2,6 @@
import sys
-import os
+
+from mercurial import ui, hg
if len(sys.argv)>1 and sys.argv[1] in ["-h","--help"]:
@@ -10,31 +11,19 @@
"""
exit(0)
-plist = os.popen("HGRCPATH='' hg parents --template='{rev}\n'").readlines()
-if len(plist)>1:
- print "You are in the process of merging"
- exit(1)
-PAR = int(plist[0])
-f = os.popen("HGRCPATH='' hg log -r 0:tip --template='{rev} {parents}\n'").\
- readlines()
-REV = -1
-lengths=[]
-for l in f:
- REV+=1
- s = l.split()
- rev = int(s[0])
- if REV != rev:
- print "Something is seriously wrong"
- exit(1)
- if len(s) == 1:
- par1 = par2 = rev - 1
- elif len(s) == 2:
- par1 = par2 = int(s[1].split(":")[0])
+u = ui.ui()
+r = hg.repository(u, ".")
+N = r.changectx(".").rev()
+lengths=[0]*(N+1)
+for i in range(N+1):
+ p=r.changectx(i).parents()
+ if p[0]:
+ p0=lengths[p[0].rev()]
else:
- par1 = int(s[1].split(":")[0])
- par2 = int(s[2].split(":")[0])
- if rev == 0:
- lengths.append(0)
+ p0=-1
+ if len(p)>1 and p[1]:
+ p1=lengths[p[1].rev()]
else:
- lengths.append(max(lengths[par1],lengths[par2])+1)
-print lengths[PAR]
+ p1=-1
+ lengths[i]=max(p0,p1)+1
+print lengths[N]
Index: scripts/unify-sources.sh
===================================================================
--- scripts/unify-sources.sh (revision 208)
+++ scripts/unify-sources.sh (revision 411)
@@ -4,7 +4,182 @@
HGROOT=`hg root`
-function update_header() {
+# file enumaration modes
+
+function all_files() {
+ hg status -a -m -c |
+ cut -d ' ' -f 2 | grep -E '(\.(cc|h|dox)$|Makefile\.am$)' |
+ while read file; do echo $HGROOT/$file; done
+}
+
+function modified_files() {
+ hg status -a -m |
+ cut -d ' ' -f 2 | grep -E '(\.(cc|h|dox)$|Makefile\.am$)' |
+ while read file; do echo $HGROOT/$file; done
+}
+
+function changed_files() {
+ {
+ if [ -n "$HG_PARENT1" ]
+ then
+ hg status --rev $HG_PARENT1:$HG_NODE -a -m
+ fi
+ if [ -n "$HG_PARENT2" ]
+ then
+ hg status --rev $HG_PARENT2:$HG_NODE -a -m
+ fi
+ } | cut -d ' ' -f 2 | grep -E '(\.(cc|h|dox)$|Makefile\.am$)' |
+ sort | uniq |
+ while read file; do echo $HGROOT/$file; done
+}
+
+function given_files() {
+ for file in $GIVEN_FILES
+ do
+ echo $file
+ done
+}
+
+# actions
+
+function update_action() {
+ if ! diff -q $1 $2 >/dev/null
+ then
+ echo -n " [$3 updated]"
+ rm $2
+ mv $1 $2
+ CHANGED=YES
+ fi
+}
+
+function update_warning() {
+ echo -n " [$2 warning]"
+ WARNED=YES
+}
+
+function update_init() {
+ echo Update source files...
+ TOTAL_FILES=0
+ CHANGED_FILES=0
+ WARNED_FILES=0
+}
+
+function update_done() {
+ echo $CHANGED_FILES out of $TOTAL_FILES files has been changed.
+ echo $WARNED_FILES out of $TOTAL_FILES files triggered warnings.
+}
+
+function update_begin() {
+ ((TOTAL_FILES++))
+ CHANGED=NO
+ WARNED=NO
+}
+
+function update_end() {
+ if [ $CHANGED == YES ]
+ then
+ ((++CHANGED_FILES))
+ fi
+ if [ $WARNED == YES ]
+ then
+ ((++WARNED_FILES))
+ fi
+}
+
+function check_action() {
+ if [ "$3" == 'tabs' ]
+ then
+ PATTERN=$(echo -e '\t')
+ elif [ "$3" == 'trailing spaces' ]
+ then
+ PATTERN='\ +$'
+ else
+ PATTERN='*'
+ fi
+
+ if ! diff -q $1 $2 >/dev/null
+ then
+ if [ "$PATTERN" == '*' ]
+ then
+ diff $1 $2 | grep '^[0-9]' | sed "s|^\(.*\)c.*$|$2:\1: check failed: $3|g" |
+ sed "s/:\([0-9]*\),\([0-9]*\):\(.*\)$/:\1:\3 (until line \2)/g"
+ else
+ grep -n -E "$PATTERN" $2 | sed "s|^\([0-9]*\):.*$|$2:\1: check failed: $3|g"
+ fi
+ FAILED=YES
+ fi
+}
+
+function check_warning() {
+ if [ "$2" == 'long lines' ]
+ then
+ grep -n -E '.{81,}' $1 | sed "s|^\([0-9]*\):.*$|$1:\1: warning: $2|g"
+ else
+ echo "$1: warning: $2"
+ fi
+ WARNED=YES
+}
+
+function check_init() {
+ echo Check source files...
+ FAILED_FILES=0
+ WARNED_FILES=0
+ TOTAL_FILES=0
+}
+
+function check_done() {
+ echo $FAILED_FILES out of $TOTAL_FILES files has been failed.
+ echo $WARNED_FILES out of $TOTAL_FILES files triggered warnings.
+
+ if [ $WARNED_FILES -gt 0 -o $FAILED_FILES -gt 0 ]
+ then
+ if [ "$WARNING" == 'INTERACTIVE' ]
+ then
+ echo -n "Are the files with errors/warnings acceptable? (yes/no) "
+ while read answer
+ do
+ if [ "$answer" == 'yes' ]
+ then
+ return 0
+ elif [ "$answer" == 'no' ]
+ then
+ return 1
+ fi
+ echo -n "Are the files with errors/warnings acceptable? (yes/no) "
+ done
+ elif [ "$WARNING" == 'WERROR' ]
+ then
+ return 1
+ fi
+ fi
+}
+
+function check_begin() {
+ ((TOTAL_FILES++))
+ FAILED=NO
+ WARNED=NO
+}
+
+function check_end() {
+ if [ $FAILED == YES ]
+ then
+ ((++FAILED_FILES))
+ fi
+ if [ $WARNED == YES ]
+ then
+ ((++WARNED_FILES))
+ fi
+}
+
+
+
+# checks
+
+function header_check() {
+ if echo $1 | grep -q -E 'Makefile\.am$'
+ then
+ return
+ fi
+
TMP_FILE=`mktemp`
- FILE_NAME=$1
(echo "/* -*- mode: C++; indent-tabs-mode: nil; -*-
@@ -26,5 +201,5 @@
*/
"
- awk 'BEGIN { pm=0; }
+ awk 'BEGIN { pm=0; }
pm==3 { print }
/\/\* / && pm==0 { pm=1;}
@@ -32,103 +207,156 @@
/\*\// && pm==1 { pm=2;}
' $1
- ) >$TMP_FILE
-
- HEADER_CH=`diff -q $TMP_FILE $FILE_NAME >/dev/null&&echo NO||echo YES`
-
- rm $FILE_NAME
- mv $TMP_FILE $FILE_NAME
-}
-
-function update_tabs() {
+ ) >$TMP_FILE
+
+ "$ACTION"_action "$TMP_FILE" "$1" header
+}
+
+function tabs_check() {
+ if echo $1 | grep -q -v -E 'Makefile\.am$'
+ then
+ OLD_PATTERN=$(echo -e '\t')
+ NEW_PATTERN=' '
+ else
+ OLD_PATTERN=' '
+ NEW_PATTERN=$(echo -e '\t')
+ fi
TMP_FILE=`mktemp`
- FILE_NAME=$1
-
- cat $1 |
- sed -e 's/\t/ /g' >$TMP_FILE
-
- TABS_CH=`diff -q $TMP_FILE $FILE_NAME >/dev/null&&echo NO||echo YES`
-
- rm $FILE_NAME
- mv $TMP_FILE $FILE_NAME
-}
-
-function remove_trailing_space() {
+ cat $1 | sed -e "s/$OLD_PATTERN/$NEW_PATTERN/g" >$TMP_FILE
+
+ "$ACTION"_action "$TMP_FILE" "$1" 'tabs'
+}
+
+function spaces_check() {
TMP_FILE=`mktemp`
- FILE_NAME=$1
-
- cat $1 |
- sed -e 's/ \+$//g' >$TMP_FILE
-
- SPACES_CH=`diff -q $TMP_FILE $FILE_NAME >/dev/null&&echo NO||echo YES`
-
- rm $FILE_NAME
- mv $TMP_FILE $FILE_NAME
-}
-
-function long_line_test() {
- cat $1 |grep -q -E '.{81,}'
-}
-
-function update_file() {
- echo -n ' update' $i ...
-
- update_header $1
- update_tabs $1
- remove_trailing_space $1
-
- CHANGED=NO;
- if [[ $HEADER_CH = YES ]];
- then
- echo -n ' [header updated]'
- CHANGED=YES;
- fi
- if [[ $TABS_CH = YES ]];
- then
- echo -n ' [tabs removed]'
- CHANGED=YES;
- fi
- if [[ $SPACES_CH = YES ]];
- then
- echo -n ' [trailing spaces removed]'
- CHANGED=YES;
- fi
- if long_line_test $1 ;
- then
- echo -n ' [LONG LINES]'
- ((LONG_LINE_FILES++))
- fi
- echo
- if [[ $CHANGED = YES ]];
- then
- ((CHANGED_FILES++))
- fi
-}
-
-CHANGED_FILES=0
-TOTAL_FILES=0
-LONG_LINE_FILES=0
-if [ $# == 0 ]; then
- echo Update all source files...
- for i in `hg manifest|grep -E '\.(cc|h|dox)$'`
+ cat $1 | sed -e 's/ \+$//g' >$TMP_FILE
+
+ "$ACTION"_action "$TMP_FILE" "$1" 'trailing spaces'
+}
+
+function long_lines_check() {
+ if cat $1 | grep -q -E '.{81,}'
+ then
+ "$ACTION"_warning $1 'long lines'
+ fi
+}
+
+# process the file
+
+function process_file() {
+ if [ "$ACTION" == 'update' ]
+ then
+ echo -n " $ACTION $1..."
+ else
+ echo " $ACTION $1..."
+ fi
+
+ CHECKING="header tabs spaces long_lines"
+
+ "$ACTION"_begin $1
+ for check in $CHECKING
do
- update_file $HGROOT/$i
- ((TOTAL_FILES++))
+ "$check"_check $1
done
- echo ' done.'
-else
- for i in $*
+ "$ACTION"_end $1
+ if [ "$ACTION" == 'update' ]
+ then
+ echo
+ fi
+}
+
+function process_all {
+ "$ACTION"_init
+ while read file
do
- update_file $i
- ((TOTAL_FILES++))
- done
+ process_file $file
+ done < <($FILES)
+ "$ACTION"_done
+}
+
+while [ $# -gt 0 ]
+do
+
+ if [ "$1" == '--help' ] || [ "$1" == '-h' ]
+ then
+ echo -n \
+"Usage:
+ $0 [OPTIONS] [files]
+Options:
+ --dry-run|-n
+ Check the files, but do not modify them.
+ --interactive|-i
+ If --dry-run is specified and the checker emits warnings,
+ then the user is asked if the warnings should be considered
+ errors.
+ --werror|-w
+ Make all warnings into errors.
+ --all|-a
+ Check all source files in the repository.
+ --modified|-m
+ Check only the modified (and new) source files. This option is
+ useful to check the modification before making a commit.
+ --changed|-c
+ Check only the changed source files compared to the parent(s) of
+ the current hg node. This option is useful as hg hook script.
+ To automatically check all your changes before making a commit,
+ add the following section to the appropriate .hg/hgrc file.
+
+ [hooks]
+ pretxncommit.checksources = scripts/unify-sources.sh -c -n -i
+
+ --help|-h
+ Print this help message.
+ files
+ The files to check/unify. If no file names are given, the modified
+ source files will be checked/unified (just like using the
+ --modified|-m option).
+"
+ exit 0
+ elif [ "$1" == '--dry-run' ] || [ "$1" == '-n' ]
+ then
+ [ -n "$ACTION" ] && echo "Conflicting action options" >&2 && exit 1
+ ACTION=check
+ elif [ "$1" == "--all" ] || [ "$1" == '-a' ]
+ then
+ [ -n "$FILES" ] && echo "Conflicting target options" >&2 && exit 1
+ FILES=all_files
+ elif [ "$1" == "--changed" ] || [ "$1" == '-c' ]
+ then
+ [ -n "$FILES" ] && echo "Conflicting target options" >&2 && exit 1
+ FILES=changed_files
+ elif [ "$1" == "--modified" ] || [ "$1" == '-m' ]
+ then
+ [ -n "$FILES" ] && echo "Conflicting target options" >&2 && exit 1
+ FILES=modified_files
+ elif [ "$1" == "--interactive" ] || [ "$1" == "-i" ]
+ then
+ [ -n "$WARNING" ] && echo "Conflicting warning options" >&2 && exit 1
+ WARNING='INTERACTIVE'
+ elif [ "$1" == "--werror" ] || [ "$1" == "-w" ]
+ then
+ [ -n "$WARNING" ] && echo "Conflicting warning options" >&2 && exit 1
+ WARNING='WERROR'
+ elif [ $(echo x$1 | cut -c 2) == '-' ]
+ then
+ echo "Invalid option $1" >&2 && exit 1
+ else
+ [ -n "$FILES" ] && echo "Invalid option $1" >&2 && exit 1
+ GIVEN_FILES=$@
+ FILES=given_files
+ break
+ fi
+
+ shift
+done
+
+if [ -z $FILES ]
+then
+ FILES=modified_files
fi
-echo $CHANGED_FILES out of $TOTAL_FILES files has been changed.
-if [[ $LONG_LINE_FILES -gt 1 ]]; then
- echo
- echo WARNING: $LONG_LINE_FILES files contains long lines!
- echo
-elif [[ $LONG_LINE_FILES -gt 0 ]]; then
- echo
- echo WARNING: a file contains long lines!
- echo
+
+if [ -z $ACTION ]
+then
+ ACTION=update
fi
+
+process_all
Index: test/CMakeLists.txt
===================================================================
--- test/CMakeLists.txt (revision 225)
+++ test/CMakeLists.txt (revision 339)
@@ -17,4 +17,5 @@
kruskal_test
maps_test
+ max_matching_test
random_test
path_test
Index: test/Makefile.am
===================================================================
--- test/Makefile.am (revision 228)
+++ test/Makefile.am (revision 404)
@@ -1,4 +1,6 @@
EXTRA_DIST += \
- test/CMakeLists.txt
+ test/CMakeLists.txt \
+ test/min_cost_flow_test.lgf \
+ test/preflow_graph.lgf
noinst_HEADERS += \
@@ -20,6 +22,9 @@
test/kruskal_test \
test/maps_test \
+ test/max_matching_test \
test/random_test \
test/path_test \
+ test/preflow_test \
+ test/suurballe_test \
test/test_tools_fail \
test/test_tools_pass \
@@ -43,5 +48,8 @@
test_kruskal_test_SOURCES = test/kruskal_test.cc
test_maps_test_SOURCES = test/maps_test.cc
+test_max_matching_test_SOURCES = test/max_matching_test.cc
test_path_test_SOURCES = test/path_test.cc
+test_preflow_test_SOURCES = test/preflow_test.cc
+test_suurballe_test_SOURCES = test/suurballe_test.cc
test_random_test_SOURCES = test/random_test.cc
test_test_tools_fail_SOURCES = test/test_tools_fail.cc
Index: test/digraph_test.cc
===================================================================
--- test/digraph_test.cc (revision 228)
+++ test/digraph_test.cc (revision 388)
@@ -20,6 +20,5 @@
#include
#include
-//#include
-//#include
+#include