Index: Makefile.am
===================================================================
--- Makefile.am (revision 321)
+++ Makefile.am (revision 363)
@@ -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 363)
@@ -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 335)
@@ -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 367)
@@ -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 337)
@@ -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 351)
@@ -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,11 @@
This group describes general \c EPS drawing methods and special
graph exporting tools.
+*/
+
+/**
+@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 335)
+++ doc/images/grid_graph.eps (revision 335)
@@ -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 343)
@@ -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 366)
@@ -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)
@@ -29,5 +29,8 @@
lemon/dim2.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 +39,10 @@
lemon/maps.h \
lemon/math.h \
+ lemon/max_matching.h \
+ lemon/nauty_reader.h \
lemon/path.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 329)
@@ -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 361)
@@ -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 361)
@@ -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 361)
@@ -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 361)
@@ -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 360)
@@ -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 361)
@@ -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/full_graph.h
===================================================================
--- lemon/full_graph.h (revision 360)
+++ lemon/full_graph.h (revision 360)
@@ -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 360)
+++ lemon/grid_graph.h (revision 360)
@@ -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 372)
+++ lemon/hypercube_graph.h (revision 372)
@@ -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 329)
@@ -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 330)
+++ lemon/max_matching.h (revision 330)
@@ -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 359)
+++ lemon/nauty_reader.h (revision 359)
@@ -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/random.h
===================================================================
--- lemon/random.h (revision 280)
+++ lemon/random.h (revision 340)
@@ -541,8 +541,4 @@
/// @{
- ///\name Initialization
- ///
- /// @{
-
/// \brief Default constructor
///
@@ -709,10 +705,4 @@
}
- /// @}
-
- ///\name Uniform distributions
- ///
- /// @{
-
/// \brief Returns a random real number from the range [0, 1)
///
@@ -771,6 +761,4 @@
return _random_bits::IntConversion::convert(core);
}
-
- /// @}
unsigned int uinteger() {
@@ -807,8 +795,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.
@@ -817,7 +804,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.
@@ -832,11 +819,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));
}
@@ -944,5 +971,4 @@
///\name Two dimensional distributions
///
-
///@{
@@ -961,5 +987,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 373)
+++ lemon/smart_graph.h (revision 375)
@@ -68,5 +68,5 @@
typedef True NodeNumTag;
- typedef True EdgeNumTag;
+ typedef True ArcNumTag;
int nodeNum() const { return nodes.size(); }
@@ -467,6 +467,6 @@
public:
- operator Edge() const {
- return _id != -1 ? edgeFromId(_id / 2) : INVALID;
+ operator Edge() const {
+ return _id != -1 ? edgeFromId(_id / 2) : INVALID;
}
@@ -483,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; }
Index: lemon/suurballe.h
===================================================================
--- lemon/suurballe.h (revision 346)
+++ lemon/suurballe.h (revision 346)
@@ -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/unify-sources.sh
===================================================================
--- scripts/unify-sources.sh (revision 208)
+++ scripts/unify-sources.sh (revision 341)
@@ -4,7 +4,185 @@
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 [ $FAILED_FILES -gt 0 ]
+ then
+ return 1
+ elif [ $WARNED_FILES -gt 0 ]
+ then
+ if [ "$WARNING" == 'INTERACTIVE' ]
+ then
+ echo -n "Are the files with 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 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 +204,5 @@
*/
"
- awk 'BEGIN { pm=0; }
+ awk 'BEGIN { pm=0; }
pm==3 { print }
/\/\* / && pm==0 { pm=1;}
@@ -32,103 +210,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 327)
@@ -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 345)
@@ -1,4 +1,5 @@
EXTRA_DIST += \
- test/CMakeLists.txt
+ test/CMakeLists.txt \
+ test/min_cost_flow_test.lgf
noinst_HEADERS += \
@@ -20,6 +21,8 @@
test/kruskal_test \
test/maps_test \
+ test/max_matching_test \
test/random_test \
test/path_test \
+ test/suurballe_test \
test/test_tools_fail \
test/test_tools_pass \
@@ -43,5 +46,7 @@
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_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 374)
+++ test/digraph_test.cc (revision 375)
@@ -20,6 +20,5 @@
#include
#include
-//#include
-//#include
+#include
#include "test_tools.h"
@@ -319,10 +318,7 @@
checkConcept, SmartDigraph>();
}
-// { // Checking FullDigraph
-// checkConcept();
-// }
-// { // Checking HyperCubeDigraph
-// checkConcept();
-// }
+ { // Checking FullDigraph
+ checkConcept();
+ }
}
@@ -375,4 +371,36 @@
check(!g.valid(g.nodeFromId(-1)), "Wrong validity check");
check(!g.valid(g.arcFromId(-1)), "Wrong validity check");
+}
+
+void checkFullDigraph(int num) {
+ typedef FullDigraph Digraph;
+ DIGRAPH_TYPEDEFS(Digraph);
+ Digraph G(num);
+
+ checkGraphNodeList(G, num);
+ checkGraphArcList(G, num * num);
+
+ for (NodeIt n(G); n != INVALID; ++n) {
+ checkGraphOutArcList(G, n, num);
+ checkGraphInArcList(G, n, num);
+ }
+
+ checkGraphConArcList(G, num * num);
+
+ checkNodeIds(G);
+ checkArcIds(G);
+ checkGraphNodeMap(G);
+ checkGraphArcMap(G);
+
+ for (int i = 0; i < G.nodeNum(); ++i) {
+ check(G.index(G(i)) == i, "Wrong index");
+ }
+
+ for (NodeIt s(G); s != INVALID; ++s) {
+ for (NodeIt t(G); t != INVALID; ++t) {
+ Arc a = G.arc(s, t);
+ check(G.source(a) == s && G.target(a) == t, "Wrong arc lookup");
+ }
+ }
}
@@ -392,4 +420,7 @@
checkDigraphValidity();
}
+ { // Checking FullDigraph
+ checkFullDigraph(8);
+ }
}
Index: test/graph_test.cc
===================================================================
--- test/graph_test.cc (revision 374)
+++ test/graph_test.cc (revision 375)
@@ -20,6 +20,7 @@
#include
#include
-// #include
-// #include
+#include
+#include
+#include
#include "test_tools.h"
@@ -262,4 +263,51 @@
}
+void checkFullGraph(int num) {
+ typedef FullGraph Graph;
+ GRAPH_TYPEDEFS(Graph);
+
+ Graph G(num);
+ checkGraphNodeList(G, num);
+ checkGraphEdgeList(G, num * (num - 1) / 2);
+
+ for (NodeIt n(G); n != INVALID; ++n) {
+ checkGraphOutArcList(G, n, num - 1);
+ checkGraphInArcList(G, n, num - 1);
+ checkGraphIncEdgeList(G, n, num - 1);
+ }
+
+ checkGraphConArcList(G, num * (num - 1));
+ checkGraphConEdgeList(G, num * (num - 1) / 2);
+
+ checkArcDirections(G);
+
+ checkNodeIds(G);
+ checkArcIds(G);
+ checkEdgeIds(G);
+ checkGraphNodeMap(G);
+ checkGraphArcMap(G);
+ checkGraphEdgeMap(G);
+
+
+ for (int i = 0; i < G.nodeNum(); ++i) {
+ check(G.index(G(i)) == i, "Wrong index");
+ }
+
+ for (NodeIt u(G); u != INVALID; ++u) {
+ for (NodeIt v(G); v != INVALID; ++v) {
+ Edge e = G.edge(u, v);
+ Arc a = G.arc(u, v);
+ if (u == v) {
+ check(e == INVALID, "Wrong edge lookup");
+ check(a == INVALID, "Wrong arc lookup");
+ } else {
+ check((G.u(e) == u && G.v(e) == v) ||
+ (G.u(e) == v && G.v(e) == u), "Wrong edge lookup");
+ check(G.source(a) == u && G.target(a) == v, "Wrong arc lookup");
+ }
+ }
+ }
+}
+
void checkConcepts() {
{ // Checking graph components
@@ -291,12 +339,13 @@
checkConcept, SmartGraph>();
}
-// { // Checking FullGraph
-// checkConcept();
-// checkGraphIterators();
-// }
-// { // Checking GridGraph
-// checkConcept();
-// checkGraphIterators();
-// }
+ { // Checking FullGraph
+ checkConcept();
+ }
+ { // Checking GridGraph
+ checkConcept();
+ }
+ { // Checking HypercubeGraph
+ checkConcept();
+ }
}
@@ -355,47 +404,130 @@
}
-// void checkGridGraph(const GridGraph& g, int w, int h) {
-// check(g.width() == w, "Wrong width");
-// check(g.height() == h, "Wrong height");
-
-// for (int i = 0; i < w; ++i) {
-// for (int j = 0; j < h; ++j) {
-// check(g.col(g(i, j)) == i, "Wrong col");
-// check(g.row(g(i, j)) == j, "Wrong row");
-// }
-// }
-
-// for (int i = 0; i < w; ++i) {
-// for (int j = 0; j < h - 1; ++j) {
-// check(g.source(g.down(g(i, j))) == g(i, j), "Wrong down");
-// check(g.target(g.down(g(i, j))) == g(i, j + 1), "Wrong down");
-// }
-// check(g.down(g(i, h - 1)) == INVALID, "Wrong down");
-// }
-
-// for (int i = 0; i < w; ++i) {
-// for (int j = 1; j < h; ++j) {
-// check(g.source(g.up(g(i, j))) == g(i, j), "Wrong up");
-// check(g.target(g.up(g(i, j))) == g(i, j - 1), "Wrong up");
-// }
-// check(g.up(g(i, 0)) == INVALID, "Wrong up");
-// }
-
-// for (int j = 0; j < h; ++j) {
-// for (int i = 0; i < w - 1; ++i) {
-// check(g.source(g.right(g(i, j))) == g(i, j), "Wrong right");
-// check(g.target(g.right(g(i, j))) == g(i + 1, j), "Wrong right");
-// }
-// check(g.right(g(w - 1, j)) == INVALID, "Wrong right");
-// }
-
-// for (int j = 0; j < h; ++j) {
-// for (int i = 1; i < w; ++i) {
-// check(g.source(g.left(g(i, j))) == g(i, j), "Wrong left");
-// check(g.target(g.left(g(i, j))) == g(i - 1, j), "Wrong left");
-// }
-// check(g.left(g(0, j)) == INVALID, "Wrong left");
-// }
-// }
+void checkGridGraph(int width, int height) {
+ typedef GridGraph Graph;
+ GRAPH_TYPEDEFS(Graph);
+ Graph G(width, height);
+
+ check(G.width() == width, "Wrong column number");
+ check(G.height() == height, "Wrong row number");
+
+ for (int i = 0; i < width; ++i) {
+ for (int j = 0; j < height; ++j) {
+ check(G.col(G(i, j)) == i, "Wrong column");
+ check(G.row(G(i, j)) == j, "Wrong row");
+ check(G.pos(G(i, j)).x == i, "Wrong column");
+ check(G.pos(G(i, j)).y == j, "Wrong row");
+ }
+ }
+
+ for (int j = 0; j < height; ++j) {
+ for (int i = 0; i < width - 1; ++i) {
+ check(G.source(G.right(G(i, j))) == G(i, j), "Wrong right");
+ check(G.target(G.right(G(i, j))) == G(i + 1, j), "Wrong right");
+ }
+ check(G.right(G(width - 1, j)) == INVALID, "Wrong right");
+ }
+
+ for (int j = 0; j < height; ++j) {
+ for (int i = 1; i < width; ++i) {
+ check(G.source(G.left(G(i, j))) == G(i, j), "Wrong left");
+ check(G.target(G.left(G(i, j))) == G(i - 1, j), "Wrong left");
+ }
+ check(G.left(G(0, j)) == INVALID, "Wrong left");
+ }
+
+ for (int i = 0; i < width; ++i) {
+ for (int j = 0; j < height - 1; ++j) {
+ check(G.source(G.up(G(i, j))) == G(i, j), "Wrong up");
+ check(G.target(G.up(G(i, j))) == G(i, j + 1), "Wrong up");
+ }
+ check(G.up(G(i, height - 1)) == INVALID, "Wrong up");
+ }
+
+ for (int i = 0; i < width; ++i) {
+ for (int j = 1; j < height; ++j) {
+ check(G.source(G.down(G(i, j))) == G(i, j), "Wrong down");
+ check(G.target(G.down(G(i, j))) == G(i, j - 1), "Wrong down");
+ }
+ check(G.down(G(i, 0)) == INVALID, "Wrong down");
+ }
+
+ checkGraphNodeList(G, width * height);
+ checkGraphEdgeList(G, width * (height - 1) + (width - 1) * height);
+ checkGraphArcList(G, 2 * (width * (height - 1) + (width - 1) * height));
+
+ for (NodeIt n(G); n != INVALID; ++n) {
+ int nb = 4;
+ if (G.col(n) == 0) --nb;
+ if (G.col(n) == width - 1) --nb;
+ if (G.row(n) == 0) --nb;
+ if (G.row(n) == height - 1) --nb;
+
+ checkGraphOutArcList(G, n, nb);
+ checkGraphInArcList(G, n, nb);
+ checkGraphIncEdgeList(G, n, nb);
+ }
+
+ checkArcDirections(G);
+
+ checkGraphConArcList(G, 2 * (width * (height - 1) + (width - 1) * height));
+ checkGraphConEdgeList(G, width * (height - 1) + (width - 1) * height);
+
+ checkNodeIds(G);
+ checkArcIds(G);
+ checkEdgeIds(G);
+ checkGraphNodeMap(G);
+ checkGraphArcMap(G);
+ checkGraphEdgeMap(G);
+
+}
+
+void checkHypercubeGraph(int dim) {
+ GRAPH_TYPEDEFS(HypercubeGraph);
+
+ HypercubeGraph G(dim);
+ checkGraphNodeList(G, 1 << dim);
+ checkGraphEdgeList(G, dim * (1 << (dim-1)));
+ checkGraphArcList(G, dim * (1 << dim));
+
+ Node n = G.nodeFromId(dim);
+
+ for (NodeIt n(G); n != INVALID; ++n) {
+ checkGraphIncEdgeList(G, n, dim);
+ for (IncEdgeIt e(G, n); e != INVALID; ++e) {
+ check( (G.u(e) == n &&
+ G.id(G.v(e)) == (G.id(n) ^ (1 << G.dimension(e)))) ||
+ (G.v(e) == n &&
+ G.id(G.u(e)) == (G.id(n) ^ (1 << G.dimension(e)))),
+ "Wrong edge or wrong dimension");
+ }
+
+ checkGraphOutArcList(G, n, dim);
+ for (OutArcIt a(G, n); a != INVALID; ++a) {
+ check(G.source(a) == n &&
+ G.id(G.target(a)) == (G.id(n) ^ (1 << G.dimension(a))),
+ "Wrong arc or wrong dimension");
+ }
+
+ checkGraphInArcList(G, n, dim);
+ for (InArcIt a(G, n); a != INVALID; ++a) {
+ check(G.target(a) == n &&
+ G.id(G.source(a)) == (G.id(n) ^ (1 << G.dimension(a))),
+ "Wrong arc or wrong dimension");
+ }
+ }
+
+ checkGraphConArcList(G, (1 << dim) * dim);
+ checkGraphConEdgeList(G, dim * (1 << (dim-1)));
+
+ checkArcDirections(G);
+
+ checkNodeIds(G);
+ checkArcIds(G);
+ checkEdgeIds(G);
+ checkGraphNodeMap(G);
+ checkGraphArcMap(G);
+ checkGraphEdgeMap(G);
+}
void checkGraphs() {
@@ -412,15 +544,21 @@
checkGraphValidity();
}
-// { // Checking FullGraph
-// FullGraph g(5);
-// checkGraphNodeList(g, 5);
-// checkGraphEdgeList(g, 10);
-// }
-// { // Checking GridGraph
-// GridGraph g(5, 6);
-// checkGraphNodeList(g, 30);
-// checkGraphEdgeList(g, 49);
-// checkGridGraph(g, 5, 6);
-// }
+ { // Checking FullGraph
+ checkFullGraph(7);
+ checkFullGraph(8);
+ }
+ { // Checking GridGraph
+ checkGridGraph(5, 8);
+ checkGridGraph(8, 5);
+ checkGridGraph(5, 5);
+ checkGridGraph(0, 0);
+ checkGridGraph(1, 1);
+ }
+ { // Checking HypercubeGraph
+ checkHypercubeGraph(1);
+ checkHypercubeGraph(2);
+ checkHypercubeGraph(3);
+ checkHypercubeGraph(4);
+ }
}
Index: test/max_matching_test.cc
===================================================================
--- test/max_matching_test.cc (revision 327)
+++ test/max_matching_test.cc (revision 327)
@@ -0,0 +1,310 @@
+/* -*- 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.
+ *
+ */
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include "test_tools.h"
+
+using namespace std;
+using namespace lemon;
+
+GRAPH_TYPEDEFS(SmartGraph);
+
+
+const int lgfn = 3;
+const std::string lgf[lgfn] = {
+ "@nodes\n"
+ "label\n"
+ "0\n"
+ "1\n"
+ "2\n"
+ "3\n"
+ "4\n"
+ "5\n"
+ "6\n"
+ "7\n"
+ "@edges\n"
+ " label weight\n"
+ "7 4 0 984\n"
+ "0 7 1 73\n"
+ "7 1 2 204\n"
+ "2 3 3 583\n"
+ "2 7 4 565\n"
+ "2 1 5 582\n"
+ "0 4 6 551\n"
+ "2 5 7 385\n"
+ "1 5 8 561\n"
+ "5 3 9 484\n"
+ "7 5 10 904\n"
+ "3 6 11 47\n"
+ "7 6 12 888\n"
+ "3 0 13 747\n"
+ "6 1 14 310\n",
+
+ "@nodes\n"
+ "label\n"
+ "0\n"
+ "1\n"
+ "2\n"
+ "3\n"
+ "4\n"
+ "5\n"
+ "6\n"
+ "7\n"
+ "@edges\n"
+ " label weight\n"
+ "2 5 0 710\n"
+ "0 5 1 241\n"
+ "2 4 2 856\n"
+ "2 6 3 762\n"
+ "4 1 4 747\n"
+ "6 1 5 962\n"
+ "4 7 6 723\n"
+ "1 7 7 661\n"
+ "2 3 8 376\n"
+ "1 0 9 416\n"
+ "6 7 10 391\n",
+
+ "@nodes\n"
+ "label\n"
+ "0\n"
+ "1\n"
+ "2\n"
+ "3\n"
+ "4\n"
+ "5\n"
+ "6\n"
+ "7\n"
+ "@edges\n"
+ " label weight\n"
+ "6 2 0 553\n"
+ "0 7 1 653\n"
+ "6 3 2 22\n"
+ "4 7 3 846\n"
+ "7 2 4 981\n"
+ "7 6 5 250\n"
+ "5 2 6 539\n",
+};
+
+void checkMatching(const SmartGraph& graph,
+ const MaxMatching& mm) {
+ int num = 0;
+
+ IntNodeMap comp_index(graph);
+ UnionFind comp(comp_index);
+
+ int barrier_num = 0;
+
+ for (NodeIt n(graph); n != INVALID; ++n) {
+ check(mm.decomposition(n) == MaxMatching::EVEN ||
+ mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
+ if (mm.decomposition(n) == MaxMatching::ODD) {
+ ++barrier_num;
+ } else {
+ comp.insert(n);
+ }
+ }
+
+ for (EdgeIt e(graph); e != INVALID; ++e) {
+ if (mm.matching(e)) {
+ check(e == mm.matching(graph.u(e)), "Wrong matching");
+ check(e == mm.matching(graph.v(e)), "Wrong matching");
+ ++num;
+ }
+ check(mm.decomposition(graph.u(e)) != MaxMatching::EVEN ||
+ mm.decomposition(graph.v(e)) != MaxMatching::MATCHED,
+ "Wrong Gallai-Edmonds decomposition");
+
+ check(mm.decomposition(graph.v(e)) != MaxMatching::EVEN ||
+ mm.decomposition(graph.u(e)) != MaxMatching::MATCHED,
+ "Wrong Gallai-Edmonds decomposition");
+
+ if (mm.decomposition(graph.u(e)) != MaxMatching::ODD &&
+ mm.decomposition(graph.v(e)) != MaxMatching::ODD) {
+ comp.join(graph.u(e), graph.v(e));
+ }
+ }
+
+ std::set comp_root;
+ int odd_comp_num = 0;
+ for (NodeIt n(graph); n != INVALID; ++n) {
+ if (mm.decomposition(n) != MaxMatching::ODD) {
+ int root = comp.find(n);
+ if (comp_root.find(root) == comp_root.end()) {
+ comp_root.insert(root);
+ if (comp.size(n) % 2 == 1) {
+ ++odd_comp_num;
+ }
+ }
+ }
+ }
+
+ check(mm.matchingSize() == num, "Wrong matching");
+ check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
+ "Wrong matching");
+ return;
+}
+
+void checkWeightedMatching(const SmartGraph& graph,
+ const SmartGraph::EdgeMap& weight,
+ const MaxWeightedMatching& mwm) {
+ for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
+ if (graph.u(e) == graph.v(e)) continue;
+ int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
+
+ for (int i = 0; i < mwm.blossomNum(); ++i) {
+ bool s = false, t = false;
+ for (MaxWeightedMatching::BlossomIt n(mwm, i);
+ n != INVALID; ++n) {
+ if (graph.u(e) == n) s = true;
+ if (graph.v(e) == n) t = true;
+ }
+ if (s == true && t == true) {
+ rw += mwm.blossomValue(i);
+ }
+ }
+ rw -= weight[e] * mwm.dualScale;
+
+ check(rw >= 0, "Negative reduced weight");
+ check(rw == 0 || !mwm.matching(e),
+ "Non-zero reduced weight on matching edge");
+ }
+
+ int pv = 0;
+ for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+ if (mwm.matching(n) != INVALID) {
+ check(mwm.nodeValue(n) >= 0, "Invalid node value");
+ pv += weight[mwm.matching(n)];
+ SmartGraph::Node o = graph.target(mwm.matching(n));
+ check(mwm.mate(n) == o, "Invalid matching");
+ check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
+ "Invalid matching");
+ } else {
+ check(mwm.mate(n) == INVALID, "Invalid matching");
+ check(mwm.nodeValue(n) == 0, "Invalid matching");
+ }
+ }
+
+ int dv = 0;
+ for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+ dv += mwm.nodeValue(n);
+ }
+
+ for (int i = 0; i < mwm.blossomNum(); ++i) {
+ check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
+ check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
+ dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
+ }
+
+ check(pv * mwm.dualScale == dv * 2, "Wrong duality");
+
+ return;
+}
+
+void checkWeightedPerfectMatching(const SmartGraph& graph,
+ const SmartGraph::EdgeMap& weight,
+ const MaxWeightedPerfectMatching& mwpm) {
+ for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
+ if (graph.u(e) == graph.v(e)) continue;
+ int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
+
+ for (int i = 0; i < mwpm.blossomNum(); ++i) {
+ bool s = false, t = false;
+ for (MaxWeightedPerfectMatching::BlossomIt n(mwpm, i);
+ n != INVALID; ++n) {
+ if (graph.u(e) == n) s = true;
+ if (graph.v(e) == n) t = true;
+ }
+ if (s == true && t == true) {
+ rw += mwpm.blossomValue(i);
+ }
+ }
+ rw -= weight[e] * mwpm.dualScale;
+
+ check(rw >= 0, "Negative reduced weight");
+ check(rw == 0 || !mwpm.matching(e),
+ "Non-zero reduced weight on matching edge");
+ }
+
+ int pv = 0;
+ for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+ check(mwpm.matching(n) != INVALID, "Non perfect");
+ pv += weight[mwpm.matching(n)];
+ SmartGraph::Node o = graph.target(mwpm.matching(n));
+ check(mwpm.mate(n) == o, "Invalid matching");
+ check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
+ "Invalid matching");
+ }
+
+ int dv = 0;
+ for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+ dv += mwpm.nodeValue(n);
+ }
+
+ for (int i = 0; i < mwpm.blossomNum(); ++i) {
+ check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
+ check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
+ dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
+ }
+
+ check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
+
+ return;
+}
+
+
+int main() {
+
+ for (int i = 0; i < lgfn; ++i) {
+ SmartGraph graph;
+ SmartGraph::EdgeMap weight(graph);
+
+ istringstream lgfs(lgf[i]);
+ graphReader(graph, lgfs).
+ edgeMap("weight", weight).run();
+
+ MaxMatching mm(graph);
+ mm.run();
+ checkMatching(graph, mm);
+
+ MaxWeightedMatching mwm(graph, weight);
+ mwm.run();
+ checkWeightedMatching(graph, weight, mwm);
+
+ MaxWeightedPerfectMatching mwpm(graph, weight);
+ bool perfect = mwpm.run();
+
+ check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
+ "Perfect matching found");
+
+ if (perfect) {
+ checkWeightedPerfectMatching(graph, weight, mwpm);
+ }
+ }
+
+ return 0;
+}
Index: test/min_cost_flow_test.lgf
===================================================================
--- test/min_cost_flow_test.lgf (revision 345)
+++ test/min_cost_flow_test.lgf (revision 345)
@@ -0,0 +1,44 @@
+@nodes
+label supply1 supply2 supply3
+1 0 20 27
+2 0 -4 0
+3 0 0 0
+4 0 0 0
+5 0 9 0
+6 0 -6 0
+7 0 0 0
+8 0 0 0
+9 0 3 0
+10 0 -2 0
+11 0 0 0
+12 0 -20 -27
+
+@arcs
+ cost capacity lower1 lower2
+1 2 70 11 0 8
+1 3 150 3 0 1
+1 4 80 15 0 2
+2 8 80 12 0 0
+3 5 140 5 0 3
+4 6 60 10 0 1
+4 7 80 2 0 0
+4 8 110 3 0 0
+5 7 60 14 0 0
+5 11 120 12 0 0
+6 3 0 3 0 0
+6 9 140 4 0 0
+6 10 90 8 0 0
+7 1 30 5 0 0
+8 12 60 16 0 4
+9 12 50 6 0 0
+10 12 70 13 0 5
+10 2 100 7 0 0
+10 7 60 10 0 0
+11 10 20 14 0 6
+12 11 30 10 0 0
+
+@attributes
+source 1
+target 12
+
+@end
Index: test/suurballe_test.cc
===================================================================
--- test/suurballe_test.cc (revision 346)
+++ test/suurballe_test.cc (revision 346)
@@ -0,0 +1,160 @@
+/* -*- 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.
+ *
+ */
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include "test_tools.h"
+
+using namespace lemon;
+
+// Check the feasibility of the flow
+template
+bool checkFlow( const Digraph& gr, const FlowMap& flow,
+ typename Digraph::Node s, typename Digraph::Node t,
+ int value )
+{
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+ for (ArcIt e(gr); e != INVALID; ++e)
+ if (!(flow[e] == 0 || flow[e] == 1)) return false;
+
+ for (NodeIt n(gr); n != INVALID; ++n) {
+ int sum = 0;
+ for (OutArcIt e(gr, n); e != INVALID; ++e)
+ sum += flow[e];
+ for (InArcIt e(gr, n); e != INVALID; ++e)
+ sum -= flow[e];
+ if (n == s && sum != value) return false;
+ if (n == t && sum != -value) return false;
+ if (n != s && n != t && sum != 0) return false;
+ }
+
+ return true;
+}
+
+// Check the optimalitiy of the flow
+template < typename Digraph, typename CostMap,
+ typename FlowMap, typename PotentialMap >
+bool checkOptimality( const Digraph& gr, const CostMap& cost,
+ const FlowMap& flow, const PotentialMap& pi )
+{
+ // Check the "Complementary Slackness" optimality condition
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+ bool opt = true;
+ for (ArcIt e(gr); e != INVALID; ++e) {
+ typename CostMap::Value red_cost =
+ cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
+ opt = (flow[e] == 0 && red_cost >= 0) ||
+ (flow[e] == 1 && red_cost <= 0);
+ if (!opt) break;
+ }
+ return opt;
+}
+
+// Check a path
+template
+bool checkPath( const Digraph& gr, const Path& path,
+ typename Digraph::Node s, typename Digraph::Node t)
+{
+ // Check the "Complementary Slackness" optimality condition
+ TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+ Node n = s;
+ for (int i = 0; i < path.length(); ++i) {
+ if (gr.source(path.nth(i)) != n) return false;
+ n = gr.target(path.nth(i));
+ }
+ return n == t;
+}
+
+
+int main()
+{
+ DIGRAPH_TYPEDEFS(ListDigraph);
+
+ // Read the test digraph
+ ListDigraph digraph;
+ ListDigraph::ArcMap length(digraph);
+ Node source, target;
+
+ std::string fname;
+ if(getenv("srcdir"))
+ fname = std::string(getenv("srcdir"));
+ else fname = ".";
+ fname += "/test/min_cost_flow_test.lgf";
+
+ std::ifstream input(fname.c_str());
+ check(input, "Input file '" << fname << "' not found");
+ DigraphReader(digraph, input).
+ arcMap("cost", length).
+ node("source", source).
+ node("target", target).
+ run();
+ input.close();
+
+ // Find 2 paths
+ {
+ Suurballe suurballe(digraph, length, source, target);
+ check(suurballe.run(2) == 2, "Wrong number of paths");
+ check(checkFlow(digraph, suurballe.flowMap(), source, target, 2),
+ "The flow is not feasible");
+ check(suurballe.totalLength() == 510, "The flow is not optimal");
+ check(checkOptimality(digraph, length, suurballe.flowMap(),
+ suurballe.potentialMap()),
+ "Wrong potentials");
+ for (int i = 0; i < suurballe.pathNum(); ++i)
+ check(checkPath(digraph, suurballe.path(i), source, target),
+ "Wrong path");
+ }
+
+ // Find 3 paths
+ {
+ Suurballe suurballe(digraph, length, source, target);
+ check(suurballe.run(3) == 3, "Wrong number of paths");
+ check(checkFlow(digraph, suurballe.flowMap(), source, target, 3),
+ "The flow is not feasible");
+ check(suurballe.totalLength() == 1040, "The flow is not optimal");
+ check(checkOptimality(digraph, length, suurballe.flowMap(),
+ suurballe.potentialMap()),
+ "Wrong potentials");
+ for (int i = 0; i < suurballe.pathNum(); ++i)
+ check(checkPath(digraph, suurballe.path(i), source, target),
+ "Wrong path");
+ }
+
+ // Find 5 paths (only 3 can be found)
+ {
+ Suurballe suurballe(digraph, length, source, target);
+ check(suurballe.run(5) == 3, "Wrong number of paths");
+ check(checkFlow(digraph, suurballe.flowMap(), source, target, 3),
+ "The flow is not feasible");
+ check(suurballe.totalLength() == 1040, "The flow is not optimal");
+ check(checkOptimality(digraph, length, suurballe.flowMap(),
+ suurballe.potentialMap()),
+ "Wrong potentials");
+ for (int i = 0; i < suurballe.pathNum(); ++i)
+ check(checkPath(digraph, suurballe.path(i), source, target),
+ "Wrong path");
+ }
+
+ return 0;
+}
Index: tools/lemon-0.x-to-1.x.sh
===================================================================
--- tools/lemon-0.x-to-1.x.sh (revision 310)
+++ tools/lemon-0.x-to-1.x.sh (revision 366)
@@ -4,124 +4,94 @@
if [ $# -eq 0 -o x$1 = "x-h" -o x$1 = "x-help" -o x$1 = "x--help" ]; then
- echo "Usage:"
- echo " $0 source-file"
- exit
+ echo "Usage:"
+ echo " $0 source-file(s)"
+ exit
fi
-TMP=`mktemp`
-
-sed -e "s/undirected graph/_gr_aph_label_/g"\
- -e "s/undirected edge/_ed_ge_label_/g"\
- -e "s/graph_/_gr_aph_label__/g"\
- -e "s/_graph/__gr_aph_label_/g"\
- -e "s/UGraph/_Gr_aph_label_/g"\
- -e "s/uGraph/_gr_aph_label_/g"\
- -e "s/ugraph/_gr_aph_label_/g"\
- -e "s/Graph/_Digr_aph_label_/g"\
- -e "s/graph/_digr_aph_label_/g"\
- -e "s/UEdge/_Ed_ge_label_/g"\
- -e "s/uEdge/_ed_ge_label_/g"\
- -e "s/uedge/_ed_ge_label_/g"\
- -e "s/IncEdgeIt/_In_cEd_geIt_label_/g"\
- -e "s/Edge/_Ar_c_label_/g"\
- -e "s/edge/_ar_c_label_/g"\
- -e "s/ANode/_Re_d_label_/g"\
- -e "s/BNode/_Blu_e_label_/g"\
- -e "s/A-Node/_Re_d_label_/g"\
- -e "s/B-Node/_Blu_e_label_/g"\
- -e "s/anode/_re_d_label_/g"\
- -e "s/bnode/_blu_e_label_/g"\
- -e "s/aNode/_re_d_label_/g"\
- -e "s/bNode/_blu_e_label_/g"\
- -e "s/_Digr_aph_label_/Digraph/g"\
- -e "s/_digr_aph_label_/digraph/g"\
- -e "s/_Gr_aph_label_/Graph/g"\
- -e "s/_gr_aph_label_/graph/g"\
- -e "s/_Ar_c_label_/Arc/g"\
- -e "s/_ar_c_label_/arc/g"\
- -e "s/_Ed_ge_label_/Edge/g"\
- -e "s/_ed_ge_label_/edge/g"\
- -e "s/_In_cEd_geIt_label_/IncEdgeIt/g"\
- -e "s/_Re_d_label_/Red/g"\
- -e "s/_Blu_e_label_/Blue/g"\
- -e "s/_re_d_label_/red/g"\
- -e "s/_blu_e_label_/blue/g"\
- -e "s/\(\W\)DefPredMap\(\W\)/\1SetPredMap\2/g"\
- -e "s/\(\W\)DefPredMap$/\1SetPredMap/g"\
- -e "s/^DefPredMap\(\W\)/SetPredMap\1/g"\
- -e "s/^DefPredMap$/SetPredMap/g"\
- -e "s/\(\W\)DefDistMap\(\W\)/\1SetDistMap\2/g"\
- -e "s/\(\W\)DefDistMap$/\1SetDistMap/g"\
- -e "s/^DefDistMap\(\W\)/SetDistMap\1/g"\
- -e "s/^DefDistMap$/SetDistMap/g"\
- -e "s/\(\W\)DefReachedMap\(\W\)/\1SetReachedMap\2/g"\
- -e "s/\(\W\)DefReachedMap$/\1SetReachedMap/g"\
- -e "s/^DefReachedMap\(\W\)/SetReachedMap\1/g"\
- -e "s/^DefReachedMap$/SetReachedMap/g"\
- -e "s/\(\W\)DefProcessedMap\(\W\)/\1SetProcessedMap\2/g"\
- -e "s/\(\W\)DefProcessedMap$/\1SetProcessedMap/g"\
- -e "s/^DefProcessedMap\(\W\)/SetProcessedMap\1/g"\
- -e "s/^DefProcessedMap$/SetProcessedMap/g"\
- -e "s/\(\W\)DefHeap\(\W\)/\1SetHeap\2/g"\
- -e "s/\(\W\)DefHeap$/\1SetHeap/g"\
- -e "s/^DefHeap\(\W\)/SetHeap\1/g"\
- -e "s/^DefHeap$/SetHeap/g"\
- -e "s/\(\W\)DefStandardHeap\(\W\)/\1SetStandradHeap\2/g"\
- -e "s/\(\W\)DefStandardHeap$/\1SetStandradHeap/g"\
- -e "s/^DefStandardHeap\(\W\)/SetStandradHeap\1/g"\
- -e "s/^DefStandardHeap$/SetStandradHeap/g"\
- -e "s/\(\W\)DefOperationTraits\(\W\)/\1SetOperationTraits\2/g"\
- -e "s/\(\W\)DefOperationTraits$/\1SetOperationTraits/g"\
- -e "s/^DefOperationTraits\(\W\)/SetOperationTraits\1/g"\
- -e "s/^DefOperationTraits$/SetOperationTraits/g"\
- -e "s/\(\W\)DefProcessedMapToBeDefaultMap\(\W\)/\1SetStandardProcessedMap\2/g"\
- -e "s/\(\W\)DefProcessedMapToBeDefaultMap$/\1SetStandardProcessedMap/g"\
- -e "s/^DefProcessedMapToBeDefaultMap\(\W\)/SetStandardProcessedMap\1/g"\
- -e "s/^DefProcessedMapToBeDefaultMap$/SetStandardProcessedMap/g"\
- -e "s/\(\W\)IntegerMap\(\W\)/\1RangeMap\2/g"\
- -e "s/\(\W\)IntegerMap$/\1RangeMap/g"\
- -e "s/^IntegerMap\(\W\)/RangeMap\1/g"\
- -e "s/^IntegerMap$/RangeMap/g"\
- -e "s/\(\W\)integerMap\(\W\)/\1rangeMap\2/g"\
- -e "s/\(\W\)integerMap$/\1rangeMap/g"\
- -e "s/^integerMap\(\W\)/rangeMap\1/g"\
- -e "s/^integerMap$/rangeMap/g"\
- -e "s/\(\W\)copyGraph\(\W\)/\1graphCopy\2/g"\
- -e "s/\(\W\)copyGraph$/\1graphCopy/g"\
- -e "s/^copyGraph\(\W\)/graphCopy\1/g"\
- -e "s/^copyGraph$/graphCopy/g"\
- -e "s/\(\W\)copyDigraph\(\W\)/\1digraphCopy\2/g"\
- -e "s/\(\W\)copyDigraph$/\1digraphCopy/g"\
- -e "s/^copyDigraph\(\W\)/digraphCopy\1/g"\
- -e "s/^copyDigraph$/digraphCopy/g"\
- -e "s/\(\W\)\([sS]\)tdMap\(\W\)/\1\2parseMap\3/g"\
- -e "s/\(\W\)\([sS]\)tdMap$/\1\2parseMap/g"\
- -e "s/^\([sS]\)tdMap\(\W\)/\1parseMap\2/g"\
- -e "s/^\([sS]\)tdMap$/\1parseMap/g"\
- -e "s/\(\W\)\([Ff]\)unctorMap\(\W\)/\1\2unctorToMap\3/g"\
- -e "s/\(\W\)\([Ff]\)unctorMap$/\1\2unctorToMap/g"\
- -e "s/^\([Ff]\)unctorMap\(\W\)/\1unctorToMap\2/g"\
- -e "s/^\([Ff]\)unctorMap$/\1unctorToMap/g"\
- -e "s/\(\W\)\([Mm]\)apFunctor\(\W\)/\1\2apToFunctor\3/g"\
- -e "s/\(\W\)\([Mm]\)apFunctor$/\1\2apToFunctor/g"\
- -e "s/^\([Mm]\)apFunctor\(\W\)/\1apToFunctor\2/g"\
- -e "s/^\([Mm]\)apFunctor$/\1apToFunctor/g"\
- -e "s/\(\W\)\([Ff]\)orkWriteMap\(\W\)/\1\2orkMap\3/g"\
- -e "s/\(\W\)\([Ff]\)orkWriteMap$/\1\2orkMap/g"\
- -e "s/^\([Ff]\)orkWriteMap\(\W\)/\1orkMap\2/g"\
- -e "s/^\([Ff]\)orkWriteMap$/\1orkMap/g"\
- -e "s/\(\W\)StoreBoolMap\(\W\)/\1LoggerBoolMap\2/g"\
- -e "s/\(\W\)StoreBoolMap$/\1LoggerBoolMap/g"\
- -e "s/^StoreBoolMap\(\W\)/LoggerBoolMap\1/g"\
- -e "s/^StoreBoolMap$/LoggerBoolMap/g"\
- -e "s/\(\W\)storeBoolMap\(\W\)/\1loggerBoolMap\2/g"\
- -e "s/\(\W\)storeBoolMap$/\1loggerBoolMap/g"\
- -e "s/^storeBoolMap\(\W\)/loggerBoolMap\1/g"\
- -e "s/^storeBoolMap$/loggerBoolMap/g"\
- -e "s/\(\W\)BoundingBox\(\W\)/\1Box\2/g"\
- -e "s/\(\W\)BoundingBox$/\1Box/g"\
- -e "s/^BoundingBox\(\W\)/Box\1/g"\
- -e "s/^BoundingBox$/Box/g"\
-<$1 > $TMP
-
-mv $TMP $1
+for i in $@
+do
+ echo Update $i...
+ TMP=`mktemp`
+ sed -e "s/\/_gr_aph_label_/g"\
+ -e "s/\/_gr_aph_label_s/g"\
+ -e "s/\/_ed_ge_label_/g"\
+ -e "s/\/_ed_ge_label_s/g"\
+ -e "s/\/_digr_aph_label_/g"\
+ -e "s/\/_digr_aph_label_s/g"\
+ -e "s/\