Merge branch 'master' of git@git.tuebingen.mpg.de:libdai
authorJoris Mooij <jorism@marvin.jorismooij.nl>
Sun, 21 Sep 2008 07:35:08 +0000 (09:35 +0200)
committerJoris Mooij <jorism@marvin.jorismooij.nl>
Sun, 21 Sep 2008 07:35:08 +0000 (09:35 +0200)
Conflicts:

include/dai/lc.h
include/dai/mr.h
src/alldai.cpp
src/lc.cpp
tests/test.cpp

43 files changed:
ChangeLog
Makefile
Makefile.win [new file with mode: 0755]
STATUS [new file with mode: 0644]
include/dai/bipgraph.h
include/dai/bp.h
include/dai/daialg.h
include/dai/enum.h
include/dai/exceptions.h [new file with mode: 0644]
include/dai/factor.h
include/dai/factorgraph.h
include/dai/hak.h
include/dai/jtree.h
include/dai/lc.h
include/dai/mf.h
include/dai/mr.h
include/dai/prob.h
include/dai/regiongraph.h
include/dai/treeep.h
include/dai/util.h
include/dai/var.h
matlab/dai.cpp
src/alldai.cpp
src/bp.cpp
src/daialg.cpp
src/exceptions.cpp [new file with mode: 0644]
src/factorgraph.cpp
src/hak.cpp
src/jtree.cpp
src/lc.cpp
src/mf.cpp
src/mr.cpp
src/properties.cpp
src/regiongraph.cpp
src/treeep.cpp
src/util.cpp
src/x2x.cpp
tests/test.cpp
tests/testall
tests/testall.bat [new file with mode: 0755]
tests/testfast.out
tests/testregression
utils/createfg.cpp

index b1e6f67..973f7a2 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,9 @@
 libDAI-0.2.2 (2008-??-??)
 -------------------------
 
+* Now compiles also under Visual C++.
+* Added Exceptions framework.
+* Added more features to utils/createfg for creating factor graphs.
 * Pervasive change of BipartiteGraph implementation (based on an idea by
   Giuseppe Passino). BipartiteGraph no longer stores the node properties 
   (former _V1 and _V2), nor does it store a dense adjacency matrix anymore, 
@@ -44,7 +47,11 @@ libDAI-0.2.2 (2008-??-??)
   normtype is now Prob::NORMPROB by default.
 * VarSet is now implemented using a std::vector<Var> instead of a
   std::set<Var>, which yields a significant speed improvement.
+* Improved MaxSpanningTreePrims algorithm (uses boost::graph)
 * Small optimization in Diffs
+* Improved ClusterGraph implementation, yielding significant speedups
+  for the JunctionTree algorithm on large factorgraphs.
+* Improved documetation
 * Interface changes:
   - VarSet::
       stateSpace() -> states()
@@ -77,7 +84,6 @@ libDAI-0.2.2 (2008-??-??)
   - Factor::max() -> Factor::maxVal()
   - toc() in util.h now returns seconds as a double
   - VarSet::operator&&
-* Added possibility to build for Windows in Makefile
 
 
 libDAI-0.2.1 (2008-05-26)
index 154debc..5740a10 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -100,7 +100,7 @@ MEXFLAGS := $(MEXFLAGS) -DSMALLMEM
 endif
 endif
 
-HEADERS = $(INC)/bipgraph.h $(INC)/diffs.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/x2x.h
+HEADERS = $(INC)/bipgraph.h $(INC)/diffs.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/x2x.h $(INC)/exceptions.h
 
 TARGETS = tests utils $(LIB)/libdai.a example testregression doc
 ifdef WITH_MATLAB
@@ -111,8 +111,8 @@ all : $(TARGETS)
 
 matlabs : matlab/dai.$(MEXEXT) matlab/dai_readfg.$(MEXEXT) matlab/dai_writefg.$(MEXEXT) matlab/dai_removeshortloops.$(MEXEXT) matlab/dai_potstrength.$(MEXEXT)
 
-$(LIB)/libdai.a : bipgraph.o daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o $(OBJECTS)
-       ar rcs $(LIB)/libdai.a bipgraph.o daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o $(OBJECTS)
+$(LIB)/libdai.a : bipgraph.o daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o exceptions.o $(OBJECTS)
+       ar rcs $(LIB)/libdai.a bipgraph.o daialg.o alldai.o clustergraph.o factorgraph.o properties.o regiongraph.o util.o weightedgraph.o x2x.o exceptions.o $(OBJECTS)
 
 tests : tests/test
 
@@ -174,6 +174,9 @@ mr.o : $(SRC)/mr.cpp $(INC)/mr.h $(HEADERS)
 properties.o : $(SRC)/properties.cpp $(HEADERS)
        $(CC) $(CCFLAGS) -c $(SRC)/properties.cpp
 
+exceptions.o: $(SRC)/exceptions.cpp $(HEADERS)
+       $(CC) $(CCFLAGS) -c $(SRC)/exceptions.cpp
+
 alldai.o : $(SRC)/alldai.cpp $(HEADERS)
        $(CC) $(CCFLAGS) -c $(SRC)/alldai.cpp
 
@@ -218,14 +221,14 @@ matlab/matlab.o : matlab/matlab.cpp matlab/matlab.h $(HEADERS)
 # UTILS
 ########
 
-utils/createfg : utils/createfg.cpp $(HEADERS) factorgraph.o weightedgraph.o util.o bipgraph.o
-       $(CC) $(CCFLAGS) -o utils/createfg utils/createfg.cpp factorgraph.o weightedgraph.o util.o bipgraph.o $(BOOSTFLAGS)
+utils/createfg : utils/createfg.cpp $(HEADERS) factorgraph.o weightedgraph.o util.o bipgraph.o exceptions.o
+       $(CC) $(CCFLAGS) -o utils/createfg utils/createfg.cpp factorgraph.o weightedgraph.o util.o bipgraph.o exceptions.o $(BOOSTFLAGS)
 
-utils/fg2dot : utils/fg2dot.cpp $(HEADERS) factorgraph.o
-       $(CC) $(CCFLAGS) -o utils/fg2dot utils/fg2dot.cpp factorgraph.o
+utils/fg2dot : utils/fg2dot.cpp $(HEADERS) factorgraph.o exceptions.o
+       $(CC) $(CCFLAGS) -o utils/fg2dot utils/fg2dot.cpp factorgraph.o exceptions.o
 
-utils/remove_short_loops : utils/remove_short_loops.cpp $(HEADERS) factorgraph.o
-       $(CC) $(CCFLAGS) -o utils/remove_short_loops utils/remove_short_loops.cpp factorgraph.o
+utils/remove_short_loops : utils/remove_short_loops.cpp $(HEADERS) factorgraph.o exceptions.o
+       $(CC) $(CCFLAGS) -o utils/remove_short_loops utils/remove_short_loops.cpp factorgraph.o exceptions.o
 
-utils/fginfo : utils/fginfo.cpp $(HEADERS) factorgraph.o bipgraph.o
-       $(CC) $(CCFLAGS) -o utils/fginfo utils/fginfo.cpp factorgraph.o bipgraph.o
+utils/fginfo : utils/fginfo.cpp $(HEADERS) factorgraph.o bipgraph.o exceptions.o
+       $(CC) $(CCFLAGS) -o utils/fginfo utils/fginfo.cpp factorgraph.o bipgraph.o exceptions.o
diff --git a/Makefile.win b/Makefile.win
new file mode 100755 (executable)
index 0000000..8edf18d
--- /dev/null
@@ -0,0 +1,245 @@
+#   Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]\r
+#   Radboud University Nijmegen, The Netherlands\r
+#   \r
+#   This file is part of libDAI.\r
+#\r
+#   libDAI is free software; you can redistribute it and/or modify\r
+#   it under the terms of the GNU General Public License as published by\r
+#   the Free Software Foundation; either version 2 of the License, or\r
+#   (at your option) any later version.\r
+#\r
+#   libDAI is distributed in the hope that it will be useful,\r
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
+#   GNU General Public License for more details.\r
+#\r
+#   You should have received a copy of the GNU General Public License\r
+#   along with libDAI; if not, write to the Free Software\r
+#   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\r
+\r
+\r
+# Enable/disable various approximate inference methods\r
+WITH_BP=true\r
+WITH_MF=true\r
+WITH_HAK=true\r
+WITH_LC=true\r
+WITH_TREEEP=true\r
+WITH_JTREE=true\r
+WITH_MR=true\r
+# Build with debug info?\r
+DEBUG=true\r
+# Build matlab interface?\r
+#WITH_MATLAB=true\r
+# New/old matlab version?\r
+NEW_MATLAB=true\r
+# Windows or linux (default)?\r
+WINDOWS=true\r
+\r
+# Directories\r
+INC=include/dai\r
+SRC=src\r
+LIB=lib\r
+\r
+# For some reason, we have to set the library path, although it is in the environment\r
+LIBS=/LIBPATH:"C:\Program Files\Microsoft Visual Studio 9.0\VC\ATLMFC\LIB" /LIBPATH:"C:\Program Files\Microsoft Visual Studio 9.0\VC\LIB" /LIBPATH:"C:\Program Files\Microsoft SDKs\Windows\v6.0A\lib"\r
+# We use the BOOST Program Options library\r
+BOOSTLIBS=/LIBPATH:C:\boost_1_36_0\stage\lib\r
+\r
+# Compile using GNU C++ Compiler\r
+CC=cl\r
+OBJEXT=obj\r
+LIBEXT=lib\r
+\r
+# Flags for the C++ compiler\r
+CCFLAGS=/I./include /IC:\boost_1_36_0 /EHsc /Ox\r
+!IFDEF DEBUG\r
+CCFLAGS=$(CCFLAGS) /Zi /DDAI_DEBUG\r
+!ELSE\r
+CCFLAGS=$(CCFLAGS)\r
+!ENDIF\r
+\r
+!IFDEF WINDOWS\r
+CCFLAGS=$(CCFLAGS) /DWINDOWS\r
+!ENDIF\r
+\r
+!IFDEF WITH_BP\r
+CCFLAGS=$(CCFLAGS) /DWITH_BP\r
+OBJECTS=$(OBJECTS) bp.obj\r
+!ENDIF\r
+!IFDEF WITH_MF\r
+CCFLAGS=$(CCFLAGS) /DWITH_MF\r
+OBJECTS=$(OBJECTS) mf.obj\r
+!ENDIF\r
+!IFDEF WITH_HAK\r
+CCFLAGS=$(CCFLAGS) /DWITH_HAK\r
+OBJECTS=$(OBJECTS) hak.obj\r
+!ENDIF\r
+!IFDEF WITH_LC\r
+CCFLAGS=$(CCFLAGS) /DWITH_LC\r
+OBJECTS=$(OBJECTS) lc.obj\r
+!ENDIF\r
+!IFDEF WITH_TREEEP\r
+CCFLAGS=$(CCFLAGS) /DWITH_TREEEP\r
+OBJECTS=$(OBJECTS) treeep.obj\r
+!ENDIF\r
+!IFDEF WITH_JTREE\r
+CCFLAGS=$(CCFLAGS) /DWITH_JTREE\r
+OBJECTS=$(OBJECTS) jtree.obj\r
+!ENDIF\r
+!IFDEF WITH_MR\r
+CCFLAGS=$(CCFLAGS) /DWITH_MR\r
+OBJECTS=$(OBJECTS) mr.obj\r
+!ENDIF\r
+\r
+!IFDEF WITH_MATLAB\r
+# Replace the following by the directory where Matlab has been installed\r
+MATLABDIR = /opt/matlab/bin\r
+# Replace the following with the extension of compiled MEX files on this platform, e.g. .mexglx for x86\r
+MEXEXT = .mexglx\r
+MEX = $(MATLABDIR)/mex\r
+MEXFLAGS = -I.\r
+!IFDEF DEBUG\r
+MEXFLAGS = $(MEXFLAGS) -g /DDAI_DEBUG\r
+!ENDIF\r
+!IFDEF NEW_MATLAB\r
+MEXFLAGS = $(MEXFLAGS) -largeArrayDims\r
+!ELSE\r
+MEXFLAGS = $(MEXFLAGS) /DSMALLMEM\r
+!ENDIF\r
+!ENDIF\r
+\r
+HEADERS = $(INC)/bipgraph.h $(INC)/diffs.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/x2x.h $(INC)/exceptions.h\r
+\r
+TARGETS = tests utils $(LIB)/libdai.lib example.exe\r
+# testregression disabled, it uses sh and awk\r
+# doc disabled, it uses doxygen\r
+!IFDEF WITH_MATLAB\r
+TARGETS = $(TARGETS) matlabs\r
+!ENDIF\r
+all : $(TARGETS)\r
+\r
+matlabs : matlab/dai.$(MEXEXT) matlab/dai_readfg.$(MEXEXT) matlab/dai_writefg.$(MEXEXT) matlab/dai_removeshortloops.$(MEXEXT) matlab/dai_potstrength.$(MEXEXT)\r
+\r
+$(LIB)/libdai.lib : bipgraph.obj daialg.obj alldai.obj clustergraph.obj factorgraph.obj properties.obj regiongraph.obj util.obj weightedgraph.obj x2x.obj exceptions.obj $(OBJECTS)\r
+       lib /out:$(LIB)/libdai.lib bipgraph.obj daialg.obj alldai.obj clustergraph.obj factorgraph.obj properties.obj regiongraph.obj util.obj weightedgraph.obj x2x.obj exceptions.obj $(OBJECTS)\r
+\r
+tests : tests/test.exe\r
+\r
+utils : utils/createfg.exe utils/fg2dot.exe utils/remove_short_loops.exe utils/fginfo.exe\r
+\r
+testregression : tests/test\r
+       echo Testing...this can take a while...\r
+       cd tests; time ./testregression; cd ..\r
+\r
+doc : $(INC)/*.h $(SRC)/*.cpp doxygen.conf\r
+       doxygen doxygen.conf\r
+\r
+clean :\r
+       del *.obj *.ilk *.pdb *.exe matlab\*.$(MEXEXT) matlab\*.obj tests\test.exe tests\*.pdb tests\*.ilk utils\*.exe utils\*.pdb utils\*.ilk $(LIB)\libdai.lib\r
+\r
+\r
+bipgraph.obj : $(SRC)/bipgraph.cpp $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/bipgraph.cpp\r
+\r
+daialg.obj : $(SRC)/daialg.cpp $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/daialg.cpp\r
+\r
+bp.obj : $(SRC)/bp.cpp $(INC)/bp.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/bp.cpp\r
+\r
+lc.obj : $(SRC)/lc.cpp $(INC)/lc.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/lc.cpp\r
+\r
+mf.obj : $(SRC)/mf.cpp $(INC)/mf.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/mf.cpp\r
+\r
+factorgraph.obj : $(SRC)/factorgraph.cpp $(INC)/factorgraph.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/factorgraph.cpp\r
+\r
+util.obj : $(SRC)/util.cpp $(INC)/util.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/util.cpp\r
+\r
+regiongraph.obj : $(SRC)/regiongraph.cpp $(INC)/regiongraph.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/regiongraph.cpp\r
+\r
+hak.obj : $(SRC)/hak.cpp $(INC)/hak.h $(HEADERS) $(INC)/regiongraph.h\r
+       $(CC) $(CCFLAGS) -c $(SRC)/hak.cpp\r
+\r
+clustergraph.obj : $(SRC)/clustergraph.cpp $(INC)/clustergraph.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/clustergraph.cpp\r
+\r
+jtree.obj : $(SRC)/jtree.cpp $(INC)/jtree.h $(HEADERS) $(INC)/weightedgraph.h $(INC)/clustergraph.h $(INC)/regiongraph.h\r
+       $(CC) $(CCFLAGS) -c $(SRC)/jtree.cpp\r
+\r
+treeep.obj : $(SRC)/treeep.cpp $(INC)/treeep.h $(HEADERS) $(INC)/weightedgraph.h $(INC)/clustergraph.h $(INC)/regiongraph.h $(INC)/jtree.h\r
+       $(CC) $(CCFLAGS) -c $(SRC)/treeep.cpp\r
+\r
+weightedgraph.obj : $(SRC)/weightedgraph.cpp $(INC)/weightedgraph.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/weightedgraph.cpp\r
+\r
+mr.obj : $(SRC)/mr.cpp $(INC)/mr.h $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/mr.cpp\r
+\r
+properties.obj : $(SRC)/properties.cpp $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/properties.cpp\r
+\r
+exceptions.obj : $(SRC)/exceptions.cpp $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/exceptions.cpp\r
+\r
+alldai.obj : $(SRC)/alldai.cpp $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/alldai.cpp\r
+\r
+x2x.obj : $(SRC)/x2x.cpp $(HEADERS)\r
+       $(CC) $(CCFLAGS) -c $(SRC)/x2x.cpp\r
+\r
+\r
+# EXAMPLE\r
+##########\r
+\r
+example.exe : example.cpp $(HEADERS) $(LIB)/libdai.lib\r
+       $(CC) $(CCFLAGS) /Feexample example.cpp /link $(LIB)/libdai.lib $(LIBS)\r
+\r
+\r
+# TESTS\r
+########\r
+\r
+tests/test.exe : tests/test.cpp $(HEADERS) $(LIB)/libdai.lib\r
+       $(CC) $(CCFLAGS) /Fetests/test tests/test.cpp /link $(LIB)/libdai.lib $(LIBS) $(BOOSTLIBS)\r
+\r
+\r
+# MATLAB INTERFACE\r
+###################\r
+\r
+matlab/dai.$(MEXEXT) : matlab/dai.cpp $(HEADERS) matlab/matlab.obj $(LIB)/libdai.lib\r
+       $(MEX) $(MEXFLAGS) -o matlab/dai matlab/dai.cpp matlab/matlab.obj $(LIB)/libdai.lib\r
+\r
+matlab/dai_readfg.$(MEXEXT) : matlab/dai_readfg.cpp $(HEADERS) factorgraph.obj matlab/matlab.obj\r
+       $(MEX) $(MEXFLAGS) -o matlab/dai_readfg matlab/dai_readfg.cpp factorgraph.obj matlab/matlab.obj\r
+\r
+matlab/dai_writefg.$(MEXEXT) : matlab/dai_writefg.cpp $(HEADERS) factorgraph.obj matlab/matlab.obj\r
+       $(MEX) $(MEXFLAGS) -o matlab/dai_writefg matlab/dai_writefg.cpp factorgraph.obj matlab/matlab.obj\r
+\r
+matlab/dai_removeshortloops.$(MEXEXT) : matlab/dai_removeshortloops.cpp $(HEADERS) factorgraph.obj matlab/matlab.obj\r
+       $(MEX) $(MEXFLAGS) -o matlab/dai_removeshortloops matlab/dai_removeshortloops.cpp factorgraph.obj matlab/matlab.obj\r
+\r
+matlab/dai_potstrength.$(MEXEXT) : matlab/dai_potstrength.cpp $(HEADERS) matlab/matlab.obj\r
+       $(MEX) $(MEXFLAGS) -o matlab/dai_potstrength matlab/dai_potstrength.cpp matlab/matlab.obj\r
+\r
+matlab/matlab.obj : matlab/matlab.cpp matlab/matlab.h $(HEADERS)\r
+       $(MEX) $(MEXFLAGS) -outdir matlab -c matlab/matlab.cpp\r
+\r
+\r
+# UTILS\r
+########\r
+\r
+utils/createfg.exe : utils/createfg.cpp $(HEADERS) factorgraph.obj weightedgraph.obj util.obj bipgraph.obj exceptions.obj\r
+       $(CC) $(CCFLAGS) /Feutils/createfg utils/createfg.cpp factorgraph.obj weightedgraph.obj util.obj bipgraph.obj exceptions.obj /link $(LIBS) $(BOOSTLIBS)\r
+\r
+utils/fg2dot.exe : utils/fg2dot.cpp $(HEADERS) factorgraph.obj exceptions.obj\r
+       $(CC) $(CCFLAGS) /Feutils/fg2dot utils/fg2dot.cpp factorgraph.obj exceptions.obj /link $(LIBS)\r
+\r
+utils/remove_short_loops.exe : utils/remove_short_loops.cpp $(HEADERS) factorgraph.obj exceptions.obj\r
+       $(CC) $(CCFLAGS) /Feutils/remove_short_loops utils/remove_short_loops.cpp factorgraph.obj exceptions.obj /link $(LIBS)\r
+\r
+utils/fginfo.exe : utils/fginfo.cpp $(HEADERS) factorgraph.obj bipgraph.obj exceptions.obj\r
+       $(CC) $(CCFLAGS) /Feutils/fginfo utils/fginfo.cpp factorgraph.obj bipgraph.obj exceptions.obj /link $(LIBS)\r
diff --git a/STATUS b/STATUS
new file mode 100644 (file)
index 0000000..d712b2c
--- /dev/null
+++ b/STATUS
@@ -0,0 +1,39 @@
+- Idea: a FactorGraph and a RegionGraph are often equipped with
+extra properties for nodes and edges. The code to initialize those
+is often quite similar; maybe this can be abstracted to classes
+like ExtFactorGraph and ExtRegionGraph (Ext stands for Extended), e.g.
+template <typename NodeProperties, typename EdgeProperties>
+class ExtFactorGraph : public FactorGraph {
+       public:
+               std::vector<NodeProperties>               nodeProps;
+               std::vector<std::vector<EdgeProperties> > edgeProps;
+       // blabla
+}
+A disadvantage of this approach may be worse cachability.
+- http://www.boost.org/development/requirements.html#Design_and_Programming
+- Would it be a good idea to cache second-order neighborhoods (delta's) in BipGraph?
+- Would it be a good idea to add the variable label -> index hashmap in FactorGraph,
+  to replace the linear searches that are performed every time for findVar()?
+  No, a better idea is to avoid calls to findVar() as much as possible.
+- Can the FactorGraph constructors be optimized further?
+
+TODO FOR RELEASE:
+- Test Visual C++ compatibility
+- Figure out which libraries are required and document in README
+  boost headers, boost::program_options library, boost::graph library,
+  boost::math library (under Windows)
+
+FILES IN SVN HEAD THAT ARE NO LONGER RELEVANT FOR GIT MASTER
+diffs.h
+index.h
+util.h
+util.cpp
+bipgraph.h
+weightedgraph.h
+clustergraph.h
+clustergraph.cpp
+varset.h
+var.h
+utils/createfg.cpp
+exceptions.h
+exceptions.cpp
index c9701ef..19eb434 100644 (file)
@@ -241,12 +241,35 @@ class BipartiteGraph {
             _nb2.push_back( nbs2new );
         }
 
-        /// Remove node of type 1 and all incident edges.
+        /// Remove node n1 of type 1 and all incident edges.
         void erase1( size_t n1 );
 
-        /// Remove node of type 2 and all incident edges.
+        /// Remove node n2 of type 2 and all incident edges.
         void erase2( size_t n2 );
 
+        /// Add edge between node n1 of type 1 and node n2 of type 2.
+        /** If check == true, only adds the edge if it does not exist already.
+         */
+        void addEdge( size_t n1, size_t n2, bool check = true ) {
+            assert( n1 < nr1() );
+            assert( n2 < nr2() );
+            bool exists = false;
+            if( check ) {
+                // Check whether the edge already exists
+                foreach( const Neighbor &nb2, nb1(n1) )
+                    if( nb2 == n2 ) {
+                        exists = true;
+                        break;
+                    }
+            }
+            if( !exists ) { // Add edge
+                Neighbor nb_1( _nb1[n1].size(), n2, _nb2[n2].size() );
+                Neighbor nb_2( nb_1.dual, n1, nb_1.iter );
+                _nb1[n1].push_back( nb_1 );
+                _nb2[n2].push_back( nb_2 );
+            }
+        }
+
         /// Calculate second-order neighbors (i.e., neighbors of neighbors) of node n1 of type 1.
         /** If include == true, include n1 itself, otherwise exclude n1.
          */
@@ -282,11 +305,11 @@ void BipartiteGraph::create( size_t nr1, size_t nr2, EdgeInputIterator begin, Ed
     _nb2.resize( nr2 );
 
     for( EdgeInputIterator e = begin; e != end; e++ ) {
-        // Each edge yields a neighbor pair
-        Neighbor nb_1( _nb1[e->first].size(), e->second, _nb2[e->second].size() );
-        Neighbor nb_2( nb_1.dual, e->first, nb_1.iter );
-        _nb1[e->first].push_back( nb_1 );
-        _nb2[e->second].push_back( nb_2 );
+#ifdef DAI_DEBUG
+        addEdge( e->first, e->second, true );
+#else
+        addEdge( e->first, e->second, false );
+#endif
     }
 }
 
index 890b72b..04cfd7b 100644 (file)
@@ -101,7 +101,7 @@ class BP : public DAIAlgFG {
         Factor belief (const Var &n) const;
         Factor belief (const VarSet &n) const;
         std::vector<Factor> beliefs() const;
-        Complex logZ() const;
+        Real logZ() const;
 
         void init( const VarSet &ns );
         void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
index fdb0ee7..61ee8dc 100644 (file)
@@ -58,7 +58,7 @@ class InfAlg {
         virtual std::vector<Factor> beliefs() const = 0;
 
         /// Get log partition sum
-        virtual Complex logZ() const = 0;
+        virtual Real logZ() const = 0;
 
         /// Clear messages and beliefs
         virtual void init() = 0;
@@ -79,38 +79,14 @@ class InfAlg {
         /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
         virtual void clamp( const Var & n, size_t i ) = 0;
 
-        /// Return all variables that interact with var(i)
-        virtual VarSet delta( size_t i ) const = 0;
-
         /// Set all factors interacting with var(i) to 1
         virtual void makeCavity( size_t i ) = 0;
 
-        /// Get index of variable n
-        virtual size_t findVar( const Var & n ) const = 0;
-
-        /// Get index of first factor involving ns
-        virtual size_t findFactor( const VarSet &ns ) const = 0;
-
-        /// Get number of variables
-        virtual size_t nrVars() const = 0;
-
-        /// Get number of factors
-        virtual size_t nrFactors() const = 0;
-
-        /// Get const reference to variable i
-        virtual const Var & var(size_t i) const = 0;
-
-        /// Get reference to variable i
-        virtual Var & var(size_t i) = 0;
-
-        /// Get const reference to factor I
-        virtual const Factor & factor( size_t I ) const = 0;
-
-        /// Get reference to factor I
-        virtual Factor & factor( size_t I ) = 0;
+        /// Get reference to underlying FactorGraph
+        virtual FactorGraph &fg() = 0;
 
-        /// Factor I has been updated
-        virtual void updatedFactor( size_t I ) = 0;
+        /// Get const reference to underlying FactorGraph
+        virtual const FactorGraph &fg() const = 0;
 
         /// Return maximum difference between beliefs in the last pass
         virtual double maxDiff() const = 0;
@@ -142,38 +118,14 @@ class DAIAlg : public InfAlg, public T {
         /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$) (using T::clamp)
         void clamp( const Var & n, size_t i ) { T::clamp( n, i ); }
 
-        /// Return all variables that interact with var(i) (using T::delta)
-        VarSet delta( size_t i ) const { return T::delta( i ); }
-
         /// Set all factors interacting with var(i) to 1 (using T::makeCavity)
         void makeCavity( size_t i ) { T::makeCavity( i ); }
 
-        /// Get index of variable n (using T::findVar)
-        size_t findVar( const Var & n ) const { return T::findVar(n); }
-
-        /// Get index of first factor involving ns (using T::findFactor)
-        size_t findFactor( const VarSet &ns ) const { return T::findFactor(ns); }
-
-        /// Get number of variables (using T::nrFactors)
-        size_t nrVars() const { return T::nrVars(); }
-
-        /// Get number of factors (using T::nrFactors)
-        size_t nrFactors() const { return T::nrFactors(); }
-
-        /// Get const reference to variable i (using T::var)
-        const Var & var( size_t i ) const { return T::var(i); }
-
-        /// Get reference to variable i (using T::var)
-        Var & var(size_t i) { return T::var(i); }
-
-        /// Get const reference to factor I (using T::factor)
-        const Factor & factor( size_t I ) const { return T::factor(I); }
-
-        /// Get reference to factor I (using T::factor)
-        Factor & factor( size_t I ) { return T::factor(I); }
+        /// Get reference to underlying FactorGraph
+        FactorGraph &fg() { return (FactorGraph &)(*this); }
 
-        /// Factor I has been updated (using T::updatedFactor)
-        void updatedFactor( size_t I ) { T::updatedFactor(I); }
+        /// Get const reference to underlying FactorGraph
+        const FactorGraph &fg() const { return (const FactorGraph &)(*this); }
 };
 
 
index f2618f7..bc1ccb8 100644 (file)
@@ -25,6 +25,7 @@
 
 #include <cstring>
 #include <iostream>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -53,7 +54,7 @@ namespace dai {
                    break;\
                }\
            if( i == sizeof(labels) / sizeof(char const *) )\
-               throw "Unknown " #x " value";\
+               DAI_THROW(UNKNOWN_ENUM_VALUE);\
         }\
 \
         operator value () const { return v; }\
@@ -99,7 +100,7 @@ namespace dai {
                    break;\
                }\
            if( i == sizeof(labels) / sizeof(char const *) )\
-               throw "Unknown " #x " value";\
+               DAI_THROW(UNKNOWN_ENUM_VALUE);\
         }\
 \
         operator value () const { return v; }\
@@ -145,7 +146,7 @@ namespace dai {
                    break;\
                }\
            if( i == sizeof(labels) / sizeof(char const *) )\
-               throw "Unknown " #x " value";\
+               DAI_THROW(UNKNOWN_ENUM_VALUE);\
         }\
 \
         operator value () const { return v; }\
@@ -191,7 +192,7 @@ namespace dai {
                    break;\
                }\
            if( i == sizeof(labels) / sizeof(char const *) )\
-               throw "Unknown " #x " value";\
+               DAI_THROW(UNKNOWN_ENUM_VALUE);\
         }\
 \
         operator value () const { return v; }\
@@ -237,7 +238,7 @@ namespace dai {
                    break;\
                }\
            if( i == sizeof(labels) / sizeof(char const *) )\
-               throw "Unknown " #x " value";\
+               DAI_THROW(UNKNOWN_ENUM_VALUE);\
         }\
 \
         operator value () const { return v; }\
diff --git a/include/dai/exceptions.h b/include/dai/exceptions.h
new file mode 100644 (file)
index 0000000..2a5e8f4
--- /dev/null
@@ -0,0 +1,66 @@
+/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
+    Radboud University Nijmegen, The Netherlands
+    
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#ifndef __defined_libdai_exceptions_h
+#define __defined_libdai_exceptions_h
+
+
+#include <exception>
+#include <stdexcept>
+#include <string>
+
+
+#define DAI_QUOTE(x) #x
+#define DAI_TOSTRING(x) DAI_QUOTE(x)
+
+#define DAI_THROW(cod) throw dai::Exception(dai::Exception::cod, std::string(__FILE__ ", line " DAI_TOSTRING(__LINE__)))
+
+
+namespace dai {
+
+
+    class Exception : public std::runtime_error {
+        public:
+            Exception(size_t code, const std::string& msg = "") : std::runtime_error(ErrorStrings[code] + " [" +  msg + "]") {}
+
+        enum {NOT_IMPLEMENTED,
+              UNKNOWN_DAI_ALGORITHM,
+              UNKNOWN_PROPERTY_TYPE,
+              MALFORMED_PROPERTY,
+              UNKNOWN_ENUM_VALUE,
+              CANNOT_READ_FILE,
+              CANNOT_WRITE_FILE,
+              INVALID_FACTORGRAPH_FILE,
+              NOT_ALL_PROPERTIES_SPECIFIED,
+              MULTIPLE_UNDO,
+              FACTORGRAPH_NOT_CONNECTED,
+              INTERNAL_ERROR,
+              NUM_ERRORS};  // NUM_ERRORS should be the last entry
+
+        private:
+            static std::string ErrorStrings[NUM_ERRORS];
+    };
+
+
+}
+
+
+#endif
index 41b5eaf..6f2e5b0 100644 (file)
@@ -36,16 +36,15 @@ namespace dai {
 
 template<typename T> class      TFactor;
 typedef TFactor<Real>           Factor;
-typedef TFactor<Complex>        CFactor;
 
 
 // predefine friends
 template<typename T> Real           dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt );
-template<typename T> Complex        KL_dist( const TFactor<T> & p, const TFactor<T> & q );
+template<typename T> Real           KL_dist( const TFactor<T> & p, const TFactor<T> & q );
 template<typename T> std::ostream&  operator<< (std::ostream& os, const TFactor<T>& P);
 
         
-// T should be castable from and to double and to complex
+// T should be castable from and to double
 template <typename T> class TFactor {
     protected:
         VarSet      _vs;
@@ -206,13 +205,6 @@ template <typename T> class TFactor {
             return l0; 
         }
 
-        CFactor clog0() const {
-            CFactor l0; 
-            l0._vs = _vs; 
-            l0._p = _p.clog0(); 
-            return l0; 
-        }
-
         T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); }
         TFactor<T> normalized( typename Prob::NormType norm ) const { 
             TFactor<T> result;
@@ -247,7 +239,7 @@ template <typename T> class TFactor {
         T totalSum() const { return _p.totalSum(); }
         T maxAbs() const { return _p.maxAbs(); }
         T maxVal() const { return _p.maxVal(); }
-        Complex entropy() const { return _p.entropy(); }
+        Real entropy() const { return _p.entropy(); }
         T strength( const Var &i, const Var &j ) const;
 
         friend Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt ) {
@@ -260,7 +252,7 @@ template <typename T> class TFactor {
                 return dist( x._p, y._p, dt );
             }
         }
-        friend Complex KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
+        friend Real KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
         template<class U> friend std::ostream& operator<< (std::ostream& os, const TFactor<U>& P);
 };
 
@@ -302,7 +294,7 @@ template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& Q) cons
 }
 
 
-template<typename T> Complex KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
+template<typename T> Real KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
     if( P._vs.empty() || Q._vs.empty() )
         return -1;
     else {
index 77f5f8a..567d9bb 100644 (file)
@@ -25,7 +25,6 @@
 
 #include <iostream>
 #include <map>
-#include <tr1/unordered_map>
 #include <dai/bipgraph.h>
 #include <dai/factor.h>
 
@@ -48,13 +47,12 @@ class FactorGraph {
 
     protected:
         std::map<size_t,Prob>  _undoProbs;
-        Prob::NormType         _normtype;
 
     public:
         /// Default constructor
-        FactorGraph() : G(), vars(), factors(), _undoProbs(), _normtype(Prob::NORMPROB) {};
+        FactorGraph() : G(), vars(), factors(), _undoProbs() {};
         /// Copy constructor
-        FactorGraph(const FactorGraph & x) : G(x.G), vars(x.vars), factors(x.factors), _undoProbs(x._undoProbs), _normtype(x._normtype) {};
+        FactorGraph(const FactorGraph & x) : G(x.G), vars(x.vars), factors(x.factors), _undoProbs(x._undoProbs) {};
         /// Construct FactorGraph from vector of Factors
         FactorGraph(const std::vector<Factor> &P);
         // Construct a FactorGraph from given factor and variable iterators
@@ -68,7 +66,6 @@ class FactorGraph {
                 vars       = x.vars;
                 factors    = x.factors;
                 _undoProbs = x._undoProbs;
-                _normtype  = x._normtype;
             }
             return *this;
         }
@@ -129,7 +126,6 @@ class FactorGraph {
         virtual void clamp( const Var & n, size_t i );
         
         bool hasNegatives() const;
-        Prob::NormType NormType() const { return _normtype; }
         
         std::vector<VarSet> Cliques() const;
 
@@ -138,8 +134,6 @@ class FactorGraph {
         virtual void undoProb( size_t I );
         void saveProb( size_t I );
 
-        virtual void updatedFactor( size_t /*I*/ ) {};
-
     private:
         /// Part of constructors (creates edges, neighbors and adjacency matrix)
         void createGraph( size_t nrEdges );
@@ -148,7 +142,7 @@ class FactorGraph {
 
 // assumes that the set of variables in [var_begin,var_end) is the union of the variables in the factors in [fact_begin, fact_end)
 template<typename FactorInputIterator, typename VarInputIterator>
-FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : G(), _undoProbs(), _normtype(Prob::NORMPROB) {
+FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : G(), _undoProbs() {
     // add factors
     size_t nrEdges = 0;
     factors.reserve( nr_fact_hint );
index a874950..7ce573b 100644 (file)
@@ -98,7 +98,7 @@ class HAK : public DAIAlgRG {
         Factor belief( const Var &n ) const;
         Factor belief( const VarSet &ns ) const;
         std::vector<Factor> beliefs() const;
-        Complex logZ () const;
+        Real logZ () const;
 
         void init( const VarSet &ns );
         void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
index 74f76b3..9fd1637 100644 (file)
@@ -84,7 +84,7 @@ class JTree : public DAIAlgRG {
         Factor belief( const Var &n ) const;
         Factor belief( const VarSet &ns ) const;
         std::vector<Factor> beliefs() const;
-        Complex logZ() const;
+        Real logZ() const;
 
         void init( const VarSet &/*ns*/ ) {}
         void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
index 57037d9..f1fab99 100644 (file)
@@ -28,6 +28,7 @@
 #include <dai/enum.h>
 #include <dai/factorgraph.h>
 #include <dai/properties.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -92,19 +93,32 @@ class LC : public DAIAlgFG {
 
         std::string identify() const;
         Factor belief (const Var &n) const { return( _beliefs[findVar(n)] ); }
-        Factor belief (const VarSet &/*ns*/) const { assert( 0 == 1 ); }
+        Factor belief (const VarSet &/*ns*/) const { 
+            DAI_THROW(NOT_IMPLEMENTED);
+            return Factor(); 
+        }
         std::vector<Factor> beliefs() const { return _beliefs; }
-        Complex logZ() const { return NAN; }
+        Real logZ() const { 
+            DAI_THROW(NOT_IMPLEMENTED);
+            return 0.0; 
+        }
         void CalcBelief (size_t i);
         const Factor &belief (size_t i) const { return _beliefs[i]; };
         const Factor &pancake (size_t i) const { return _pancakes[i]; };
         const Factor &cavitydist (size_t i) const { return _cavitydists[i]; };
 
-        void clamp( const Var &/*n*/, size_t /*i*/ ) { assert( 0 == 1 ); }
-        void undoProbs( const VarSet &/*ns*/ ) { assert( 0 == 1 ); }
-        void saveProbs( const VarSet &/*ns*/ ) { assert( 0 == 1 ); }
-        virtual void makeCavity(const Var & /*n*/) { assert( 0 == 1 ); }
-
+        void clamp( const Var &/*n*/, size_t /*i*/ ) { 
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
+        void undoProbs( const VarSet &/*ns*/ ) { 
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
+        void saveProbs( const VarSet &/*ns*/ ) { 
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
+        virtual void makeCavity(const Var & /*n*/) { 
+            DAI_THROW(NOT_IMPLEMENTED);
+        }
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         double maxDiff() const { return maxdiff; }
index 01d6721..2befe07 100644 (file)
@@ -75,7 +75,7 @@ class MF : public DAIAlgFG {
         Factor belief (const Var &n) const;
         Factor belief (const VarSet &ns) const;
         std::vector<Factor> beliefs() const;
-        Complex logZ() const;
+        Real logZ() const;
 
         void init( const VarSet &ns );
         void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
index f3cee73..41610e7 100644 (file)
@@ -29,6 +29,7 @@
 #include <dai/daialg.h>
 #include <dai/enum.h>
 #include <dai/properties.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -77,9 +78,15 @@ class MR : public DAIAlgFG {
         void solveM();
         double run();
         Factor belief( const Var &n ) const;
-        Factor belief( const VarSet &/*ns*/ ) const { assert( 0 == 1 ); }
+        Factor belief( const VarSet &/*ns*/ ) const { 
+            DAI_THROW(NOT_IMPLEMENTED);
+            return Factor(); 
+        }
         std::vector<Factor> beliefs() const;
-        Complex logZ() const { return NAN; }
+        Real logZ() const { 
+            DAI_THROW(NOT_IMPLEMENTED);
+            return 0.0; 
+        }
         void init() {}
         static const char *Name;
         std::string identify() const;
@@ -95,7 +102,7 @@ class MR : public DAIAlgFG {
         void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
 
         double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
-        MR* clone() const { assert( 0 == 1 ); }
+        MR* clone() const { return new MR(*this); }
 
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
index 8b5f6ef..6037339 100644 (file)
@@ -23,7 +23,6 @@
 #define __defined_libdai_prob_h
 
 
-#include <complex>
 #include <cmath>
 #include <vector>
 #include <ostream>
@@ -38,15 +37,13 @@ namespace dai {
 
 
 typedef double                  Real;
-typedef std::complex<double>    Complex;
 
 template<typename T> class      TProb;
 typedef TProb<Real>             Prob;
-typedef TProb<Complex>          CProb;
 
 
 /// TProb<T> implements a probability vector of type T.
-/// T should be castable from and to double and to complex.
+/// T should be castable from and to double.
 template <typename T> class TProb {
     protected:
         /// The entries
@@ -54,7 +51,7 @@ template <typename T> class TProb {
 
     private:
         /// Calculate x times log(x), or 0 if x == 0
-        Complex xlogx( Real x ) const { return( x == 0.0 ? 0.0 : Complex(x) * std::log(Complex(x))); }
+        Real xlogx( Real x ) const { return( x == 0.0 ? 0.0 : x * std::log(x)); }
 
     public:
         /// NORMPROB means that the sum of all entries should be 1
@@ -233,9 +230,6 @@ template <typename T> class TProb {
             assert( size() == q.size() );
 #endif
             for( size_t i = 0; i < size(); i++ ) {
-#ifdef DAI_DEBUG
-//              assert( q[i] != 0.0 );
-#endif
                 if( q[i] == 0.0 )
                     _p[i] = 0.0;
                 else
@@ -334,15 +328,6 @@ template <typename T> class TProb {
             return l0;
         }
 
-        /// Pointwise (complex) log (or 0 if == 0)
-/*      CProb clog0() const {
-            CProb l0;
-            l0._p.reserve( size() );
-            for( size_t i = 0; i < size(); i++ )
-                l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( Complex( _p[i] ) ) );
-            return l0;
-        }*/
-
         /// Return distance of p and q
         friend Real dist( const TProb<T> & p, const TProb<T> & q, DistType dt ) {
 #ifdef DAI_DEBUG
@@ -372,16 +357,16 @@ template <typename T> class TProb {
             return result;
         }
 
-        /// Return (complex) Kullback-Leibler distance with q
-        friend Complex KL_dist( const TProb<T> & p, const TProb<T> & q ) {
+        /// Return Kullback-Leibler distance with q
+        friend Real KL_dist( const TProb<T> & p, const TProb<T> & q ) {
 #ifdef DAI_DEBUG
             assert( p.size() == q.size() );
 #endif
-            Complex result = 0.0;
+            Real result = 0.0;
             for( size_t i = 0; i < p.size(); i++ ) {
                 if( (Real) p[i] != 0.0 ) {
-                    Complex p_i = p[i];
-                    Complex q_i = q[i];
+                    Real p_i = p[i];
+                    Real q_i = q[i];
                     result += p_i * (std::log(p_i) - std::log(q_i));
                 }
             }
@@ -447,9 +432,9 @@ template <typename T> class TProb {
             return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less_equal<Real>(), 0.0 ) ) != _p.end());
         }
 
-        /// Returns (complex) entropy
-        Complex entropy() const {
-            Complex S = 0.0;
+        /// Returns entropy
+        Real entropy() const {
+            Real S = 0.0;
             for( size_t i = 0; i < size(); i++ )
                 S -= xlogx(_p[i]);
             return S;
index 09e17a5..16916c0 100644 (file)
@@ -189,9 +189,6 @@ class RegionGraph : public FactorGraph {
         /// We have to overload FactorGraph::undoProb because the corresponding outer regions have to be recomputed
         void undoProb( size_t I ) { FactorGraph::undoProb( I ); RecomputeOR( I ); }
 
-        /// If updateFactor is called, we know that factor I has been changed and we should recompute the outer regions involving I
-        void updatedFactor( size_t I ) { RecomputeOR( I ); }
-
         /// Send RegionGraph to output stream
         friend std::ostream & operator << ( std::ostream & os, const RegionGraph & rg );
 };
index 9b18a4c..f6a4813 100644 (file)
@@ -123,7 +123,7 @@ class TreeEP : public JTree {
         std::string identify() const;
         void init();
         double run();
-        Complex logZ() const;
+        Real logZ() const;
 
         bool offtree( size_t I ) const { return (fac2OR[I] == -1U); }
 
index d83c204..900f275 100644 (file)
 #include <cstdio>
 #include <boost/foreach.hpp>
 
+
 #ifdef WINDOWS
-    #include <float.h>  // for _isnan (?)
-    #include <boost/math/special_functions/atanh.hpp>  // for atanh
-    #include <boost/math/special_functions/log1p.hpp>  // for log1p
-//    #include <xmath.h>
+    #include <map>
+#else
+    #include <tr1/unordered_map>
 #endif
 
 
 #define foreach BOOST_FOREACH
 
 
+#ifdef WINDOWS
+    bool isnan( double x );
+    double atanh( double x );
+    double log1p( double x );
+#endif
+
+
 namespace dai {
 
 
 #ifdef WINDOWS
-    bool isnan( double x ) inline { return _isnan( x ); }
-    using boost::math:atanh;
-    using boost::math::log1p;
+    template <typename T, typename U>
+        class hash_map : public std::map<T,U> {};
+#else
+    template <typename T, typename U>
+        class hash_map : public std::tr1::unordered_map<T,U> {};
 #endif
 
 
index e9deda1..7deaa95 100644 (file)
@@ -31,7 +31,7 @@ namespace dai {
 
 
 /// Represents a discrete variable.
-/*  It contains the label of the variable, an integer-valued
+/**  It contains the label of the variable, an integer-valued
  *  unique ID for that variable, and the number of possible 
  *  values ("states") of the variable.
  */
index c4e766b..7ad11b8 100644 (file)
@@ -103,7 +103,7 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
 
 
     // Save logZ
-    double logZ = real( obj->logZ() );
+    double logZ = obj->logZ();
 
     // Save maxdiff
     double maxdiff = obj->MaxDiff();
index cca6ec2..05359ba 100644 (file)
@@ -22,6 +22,7 @@
 #include <string>
 #include <dai/alldai.h>
 #include <dai/properties.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -59,7 +60,7 @@ InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const PropertySet
     if( name == MR::Name )
         return new MR (fg, opts);
 #endif
-    throw "Unknown inference algorithm";
+    DAI_THROW(UNKNOWN_DAI_ALGORITHM);
 }
 
 
index 94680d6..0352565 100644 (file)
@@ -172,7 +172,7 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
     const ind_t ind = index(i,_I);
     for( size_t r = 0; r < prod.size(); ++r )
         marg[ind[r]] += prod[r];
-    marg.normalize( _normtype );
+    marg.normalize( Prob::NORMPROB );
     
     // Store result
     if( props.logdomain )
@@ -393,10 +393,10 @@ Factor BP::beliefF (size_t I) const {
 }
 
 
-Complex BP::logZ() const {
-    Complex sum = 0.0;
+Real BP::logZ() const {
+    Real sum = 0.0;
     for(size_t i = 0; i < nrVars(); ++i )
-        sum += Complex(1.0 - nbV(i).size()) * beliefV(i).entropy();
+        sum += (1.0 - nbV(i).size()) * beliefV(i).entropy();
     for( size_t I = 0; I < nrFactors(); ++I )
         sum -= KL_dist( beliefF(I), factor(I) );
     return sum;
index 7453f06..e199670 100644 (file)
@@ -38,7 +38,7 @@ Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) {
     if( !reInit )
         clamped->init();
 
-    Complex logZ0;
+    Real logZ0 = 0.0;
     for( State s(ns); s.valid(); s++ ) {
         // save unclamped factors connected to ns
         clamped->saveProbs( ns );
@@ -52,18 +52,16 @@ Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) {
             clamped->init();
         clamped->run();
 
-        Complex Z;
+        Real Z;
         if( s == 0 ) {
             logZ0 = clamped->logZ();
             Z = 1.0;
         } else {
             // subtract logZ0 to avoid very large numbers
             Z = exp(clamped->logZ() - logZ0);
-            if( fabs(imag(Z)) > 1e-5 )
-                cout << "Marginal:: WARNING: complex Z (" << Z << ")" << endl;
         }
 
-        Pns[s] = real(Z);
+        Pns[s] = Z;
         
         // restore clamped factors
         clamped->undoProbs( ns );
@@ -96,7 +94,7 @@ vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reIni
     if( !reInit )
         clamped->init();
 
-    Complex logZ0;
+    Real logZ0 = 0.0;
     for( size_t j = 0; j < N; j++ ) {
         // clamp Var j to its possible values
         for( size_t j_val = 0; j_val < vns[j].states(); j_val++ ) {
@@ -115,10 +113,7 @@ vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reIni
                 logZ0 = clamped->logZ();
             } else {
                 // subtract logZ0 to avoid very large numbers
-                Complex Z = exp(clamped->logZ() - logZ0);
-                if( fabs(imag(Z)) > 1e-5 )
-                    cout << "calcPairBelief::  Warning: complex Z: " << Z << endl;
-                Z_xj = real(Z);
+                Z_xj = exp(clamped->logZ() - logZ0);
             }
 
             for( size_t k = 0; k < N; k++ ) 
@@ -172,7 +167,7 @@ vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool re
     if( !reInit )
         clamped->init();
 
-    Complex logZ0;
+    Real logZ0 = 0.0;
     VarSet::const_iterator nj = ns.begin();
     for( long j = 0; j < (long)ns.size() - 1; j++, nj++ ) {
         size_t k = 0;
@@ -196,10 +191,7 @@ vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool re
                         logZ0 = clamped->logZ();
                     } else {
                         // subtract logZ0 to avoid very large numbers
-                        Complex Z = exp(clamped->logZ() - logZ0);
-                        if( fabs(imag(Z)) > 1e-5 )
-                            cout << "calcPairBelief::  Warning: complex Z: " << Z << endl;
-                        Z_xj = real(Z);
+                        Z_xj = exp(clamped->logZ() - logZ0);
                     }
 
                     // we assume that j.label() < k.label()
diff --git a/src/exceptions.cpp b/src/exceptions.cpp
new file mode 100644 (file)
index 0000000..6f8f787
--- /dev/null
@@ -0,0 +1,44 @@
+/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
+    Radboud University Nijmegen, The Netherlands
+    
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#include <dai/exceptions.h>
+
+
+namespace dai {
+
+
+    std::string Exception::ErrorStrings[NUM_ERRORS] = {
+        "This feature is not implemented",
+        "Unknown DAI algorithm",
+        "Unknown Property type",
+        "Malformed Property",
+        "Unknown ENUM value",
+        "Cannot read file",
+        "Cannot write file",
+        "Invalid FactorGraph file",
+        "Not all mandatory Properties specified",
+        "Multiple undo levels unsupported",
+        "FactorGraph is not connected",
+        "Internal error"
+    }; 
+
+
+}
index b887316..4df27c8 100644 (file)
@@ -27,9 +27,9 @@
 #include <string>
 #include <algorithm>
 #include <functional>
-#include <tr1/unordered_map>
 #include <dai/factorgraph.h>
 #include <dai/util.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -38,7 +38,7 @@ namespace dai {
 using namespace std;
 
 
-FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _undoProbs(), _normtype(Prob::NORMPROB) {
+FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _undoProbs() {
     // add factors, obtain variables
     set<Var> _vars;
     factors.reserve( P.size() );
@@ -62,7 +62,7 @@ FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _undoProbs(), _n
 /// Part of constructors (creates edges, neighbours and adjacency matrix)
 void FactorGraph::createGraph( size_t nrEdges ) {
     // create a mapping for indices
-    std::tr1::unordered_map<size_t, size_t> hashmap;
+    hash_map<size_t, size_t> hashmap;
     
     for( size_t i = 0; i < vars.size(); i++ )
         hashmap[var(i).label()] = i;
@@ -122,13 +122,13 @@ istream& operator >> (istream& is, FactorGraph& fg) {
             getline(is,line);
         is >> nr_f;
         if( is.fail() )
-            throw "ReadFromFile: unable to read number of Factors";
+            DAI_THROW(INVALID_FACTORGRAPH_FILE);
         if( verbose >= 2 )
             cout << "Reading " << nr_f << " factors..." << endl;
 
         getline (is,line);
         if( is.fail() )
-            throw "ReadFromFile: empty line expected";
+            DAI_THROW(INVALID_FACTORGRAPH_FILE);
 
         for( size_t I = 0; I < nr_f; I++ ) {
             if( verbose >= 3 )
index b4b6f0b..49903b8 100644 (file)
@@ -23,6 +23,7 @@
 #include <dai/hak.h>
 #include <dai/util.h>
 #include <dai/diffs.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -141,7 +142,7 @@ HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), props(),
                 cout << *cli << endl;
         }
     } else
-        throw "Invalid Clusters type";
+        DAI_THROW(INTERNAL_ERROR);
 
     RegionGraph rg(fg,cl);
     RegionGraph::operator=(rg);
@@ -229,10 +230,10 @@ double HAK::doGBP() {
                 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
             }
 
-            Qb_new.normalize( _normtype );
+            Qb_new.normalize( Prob::NORMPROB );
             if( Qb_new.hasNaNs() ) {
                 cout << "HAK::doGBP:  Qb_new has NaNs!" << endl;
-                return NAN;
+                return 1.0;
             }
 //          _Qb[beta] = Qb_new.makeZero(1e-100);    // damping?
             _Qb[beta] = Qb_new;
@@ -246,10 +247,10 @@ double HAK::doGBP() {
                 foreach( const Neighbor &gamma, nbOR(alpha) )
                     Qa_new *= muba(alpha,gamma.iter);
                 Qa_new ^= (1.0 / OR(alpha).c());
-                Qa_new.normalize( _normtype );
+                Qa_new.normalize( Prob::NORMPROB );
                 if( Qa_new.hasNaNs() ) {
                     cout << "HAK::doGBP:  Qa_new has NaNs!" << endl;
-                    return NAN;
+                    return 1.0;
                 }
 //              _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
                 _Qa[alpha] = Qa_new;
@@ -334,7 +335,7 @@ double HAK::doDoubleLoop() {
 
         // Inner loop
         if( isnan( doGBP() ) )
-            return NAN;
+            return 1.0;
 
         // Calculate new single variable beliefs and compare with old ones
         for( size_t i = 0; i < nrVars(); ++i ) {
@@ -419,12 +420,12 @@ vector<Factor> HAK::beliefs() const {
 }
 
 
-Complex HAK::logZ() const {
-    Complex sum = 0.0;
+Real HAK::logZ() const {
+    Real sum = 0.0;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum += Complex(IR(beta).c()) * Qb(beta).entropy();
+        sum += IR(beta).c() * Qb(beta).entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        sum += Complex(OR(alpha).c()) * Qa(alpha).entropy();
+        sum += OR(alpha).c() * Qa(alpha).entropy();
         sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
     }
     return sum;
index 39c7cfb..23be664 100644 (file)
@@ -296,12 +296,12 @@ double JTree::run() {
 }
 
 
-Complex JTree::logZ() const {
-    Complex sum = 0.0;
+Real JTree::logZ() const {
+    Real sum = 0.0;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum += Complex(IR(beta).c()) * _Qb[beta].entropy();
+        sum += IR(beta).c() * _Qb[beta].entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        sum += Complex(OR(alpha).c()) * _Qa[alpha].entropy();
+        sum += OR(alpha).c() * _Qa[alpha].entropy();
         sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
     }
     return sum;
index 2e795d4..d8da46c 100644 (file)
@@ -125,15 +125,15 @@ double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet
         cav->makeCavity( i );
 
         if( props.cavity == Properties::CavityType::FULL )
-            Bi = calcMarginal( *cav, cav->delta(i), props.reinit );
+            Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
         else if( props.cavity == Properties::CavityType::PAIR )
-            Bi = calcMarginal2ndO( *cav, cav->delta(i), props.reinit );
+            Bi = calcMarginal2ndO( *cav, cav->fg().delta(i), props.reinit );
         else if( props.cavity == Properties::CavityType::PAIR2 ) {
-            vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->delta(i), props.reinit );
+            vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->fg().delta(i), props.reinit );
             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
                 Bi *= pairbeliefs[ij];
         } else if( props.cavity == Properties::CavityType::PAIRINT ) {
-            Bi = calcMarginal( *cav, cav->delta(i), props.reinit );
+            Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
             
             // Set interactions of order > 2 to zero
             size_t N = delta(i).size();
@@ -145,7 +145,7 @@ double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet
 //            x2x::logpnorm (N, p);
             x2x::logp2p (N, p);
         } else if( props.cavity == Properties::CavityType::PAIRCUM ) {
-            Bi = calcMarginal( *cav, cav->delta(i), props.reinit );
+            Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
             
             // Set cumulants of order > 2 to zero
             size_t N = delta(i).size();
@@ -159,7 +159,7 @@ double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet
         maxdiff = cav->maxDiff();
         delete cav;
     }
-    Bi.normalize( _normtype );
+    Bi.normalize( Prob::NORMPROB );
     _cavitydists[i] = Bi;
 
     return maxdiff;
@@ -229,7 +229,7 @@ void LC::init() {
               _pancakes[i] *= _phis[i][I.iter];
         }
         
-        _pancakes[i].normalize( _normtype );
+        _pancakes[i].normalize( Prob::NORMPROB );
 
         CalcBelief(i);
     }
@@ -251,10 +251,10 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).part_sum( Ivars / var(i) );
     Factor quot = A_I.divided_by(A_Ii);
 
-    piet *= quot.divided_by( _phis[i][_I] ).normalized( _normtype );
-    _phis[i][_I] = quot.normalized( _normtype );
+    piet *= quot.divided_by( _phis[i][_I] ).normalized( Prob::NORMPROB );
+    _phis[i][_I] = quot.normalized( Prob::NORMPROB );
 
-    piet.normalize( _normtype );
+    piet.normalize( Prob::NORMPROB );
 
     if( piet.hasNaNs() ) {
         cout << "LC::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
@@ -290,7 +290,7 @@ double LC::run() {
         }
     if( hasNaNs ) {
         cout << "LC::run:  initial _pancakes has NaNs!" << endl;
-        return NAN;
+        return -1.0;
     }
 
     size_t nredges = nrEdges();
@@ -314,7 +314,7 @@ double LC::run() {
             size_t _I = update_seq[t].second;
             _pancakes[i] = NewPancake( i, _I, hasNaNs);
             if( hasNaNs )
-                return NAN;
+                return -1.0;
             CalcBelief( i );
         }
 
index 6da2bd4..efb4308 100644 (file)
@@ -109,11 +109,11 @@ double MF::run() {
             jan *= piet; 
         }
 
-        jan.normalize( _normtype );
+        jan.normalize( Prob::NORMPROB );
 
         if( jan.hasNaNs() ) {
             cout << "MF::run():  ERROR: jan has NaNs!" << endl;
-            return NAN;
+            return 1.0;
         }
 
         diffs.push( dist( jan, _beliefs[i], Prob::DISTLINF ) );
@@ -171,8 +171,8 @@ vector<Factor> MF::beliefs() const {
 }
 
 
-Complex MF::logZ() const {
-    Complex sum = 0.0;
+Real MF::logZ() const {
+    Real sum = 0.0;
     
     for(size_t i=0; i < nrVars(); i++ )
         sum -= beliefV(i).entropy();
@@ -184,7 +184,7 @@ Complex MF::logZ() const {
         Factor piet;
         piet = factor(I).log0();
         piet *= henk;
-        sum -= Complex( piet.totalSum() );
+        sum -= piet.totalSum();
     }
 
     return -sum;
index eb9ab83..b49b9c4 100644 (file)
@@ -569,7 +569,7 @@ double MR::run() {
 
         return 0.0;
     } else
-        return NAN;
+        return -1.0;
 }
 
 
index 2b131b3..d62d6d5 100644 (file)
@@ -22,6 +22,7 @@
 #include <iostream>
 #include <dai/properties.h>
 #include <dai/alldai.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -69,7 +70,7 @@ std::ostream& operator<< (std::ostream & os, const Property & p) {
         os << boost::any_cast<LC::Properties::UpdateType>(p.second);
 #endif
     else
-        throw "Unknown property type";
+        DAI_THROW(UNKNOWN_PROPERTY_TYPE);
     return( os );
 }
 
@@ -96,7 +97,7 @@ std::istream& operator >> (std::istream& is, PropertySet & ps) {
 
     // Check whether s is of the form "[.*]"
     if( (s.length() < 2) || (s.at(0) != '[') || (s.at(s.length()-1)) != ']' )
-        throw "Malformed property";
+        DAI_THROW(MALFORMED_PROPERTY);
 
     size_t N = s.length() - 1;
     for( size_t token_start = 1; token_start < N; ) {
@@ -107,7 +108,7 @@ std::istream& operator >> (std::istream& is, PropertySet & ps) {
             if( s[token_end] == '=' )
                 break;
         if( token_end == N )
-            throw "Malformed property key";
+            DAI_THROW(MALFORMED_PROPERTY);
         // we found a key
         std::string key = s.substr(token_start, token_end - token_start);
 
@@ -123,7 +124,7 @@ std::istream& operator >> (std::istream& is, PropertySet & ps) {
                 break;
         }
         if( !(level == 0) )
-            throw "Malformed property value";
+            DAI_THROW(MALFORMED_PROPERTY);
         // we found a vlue
         std::string value = s.substr(token_start, token_end - token_start);
 
index 79fb5f9..8078279 100644 (file)
@@ -114,7 +114,7 @@ RegionGraph::RegionGraph( const FactorGraph &fg, const std::vector<VarSet> &cl )
     // Create inner regions - store them in the bipartite graph
     IRs.reserve( betas.size() );
     for( set<VarSet>::const_iterator beta = betas.begin(); beta != betas.end(); beta++ )
-        IRs.push_back( Region(*beta,NAN) );
+        IRs.push_back( Region(*beta,0.0) );
     
     // Create edges
     vector<pair<size_t,size_t> > edges;
@@ -142,8 +142,9 @@ void RegionGraph::Calc_Counting_Numbers() {
     // Calculates counting numbers of inner regions based upon counting numbers of outer regions
     
     vector<vector<size_t> > ancestors(nrIRs());
+    vector<bool> assigned(nrIRs(), false);
     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
-        IR(beta).c() = NAN;
+        IR(beta).c() = 0.0;
         for( size_t beta2 = 0; beta2 < nrIRs(); beta2++ )
             if( (beta2 != beta) && IR(beta2) >> IR(beta) )
                 ancestors[beta].push_back(beta2);
@@ -153,18 +154,19 @@ void RegionGraph::Calc_Counting_Numbers() {
     do {
         new_counting = false;
         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
-            if( isnan( IR(beta).c() ) ) {
-                bool has_nan_ancestor = false;
-                for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); (beta2 != ancestors[beta].end()) && !has_nan_ancestor; beta2++ )
-                    if( isnan( IR(*beta2).c() ) )
-                        has_nan_ancestor = true;
-                if( !has_nan_ancestor ) {
+            if( !assigned[beta] ) {
+                bool has_unassigned_ancestor = false;
+                for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); (beta2 != ancestors[beta].end()) && !has_unassigned_ancestor; beta2++ )
+                    if( !assigned[*beta2] )
+                        has_unassigned_ancestor = true;
+                if( !has_unassigned_ancestor ) {
                     double c = 1.0;
                     foreach( const Neighbor &alpha, nbIR(beta) )
                         c -= OR(alpha).c();
                     for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); beta2 != ancestors[beta].end(); beta2++ )
                         c -= IR(*beta2).c();
                     IR(beta).c() = c;
+                    assigned[beta] = true;
                     new_counting = true;
                 }
             }
index 9ce7354..6f5b797 100644 (file)
@@ -214,7 +214,7 @@ TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opt
                         if( piet.vars() >> (v_i | *j) ) {
                             piet = piet.marginal( v_i | *j );
                             Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
-                            wg[UEdge(i,findVar(*j))] = real( KL_dist( piet, pietf ) );
+                            wg[UEdge(i,findVar(*j))] = KL_dist( piet, pietf );
                         } else
                             wg[UEdge(i,findVar(*j))] = 0;
                     }
@@ -251,7 +251,7 @@ TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opt
             // find maximal spanning tree
             ConstructRG( MaxSpanningTreePrims( wg ) );
         } else {
-            assert( 0 == 1 );
+            DAI_THROW(INTERNAL_ERROR);
         }
     }
 }
@@ -458,14 +458,14 @@ double TreeEP::run() {
 }
 
 
-Complex TreeEP::logZ() const {
+Real TreeEP::logZ() const {
     double sum = 0.0;
 
     // entropy of the tree
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum -= real(_Qb[beta].entropy());
+        sum -= _Qb[beta].entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        sum += real(_Qa[alpha].entropy());
+        sum += _Qa[alpha].entropy();
 
     // energy of the on-tree factors
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
index 18596ee..728f62f 100644 (file)
@@ -27,6 +27,9 @@
 
 #ifdef WINDOWS
     #include <windows.h>
+    #include <boost/math/special_functions/atanh.hpp>  // for atanh
+    #include <boost/math/special_functions/log1p.hpp>  // for log1p
+    #include <float.h>  // for _isnan
 #else
     // Assume POSIX compliant system. We need the following for querying the CPU time for this process
     #include <sys/times.h>
 #endif
 
 
+#ifdef WINDOWS
+bool isnan( double x ) {
+    return _isnan( x );
+}
+double atanh( double x ) {
+    return boost::math::atanh( x );
+}
+double log1p( double x ) {
+    return boost::math::log1p( x );
+}
+#endif
+
+
 namespace dai {
 
 
@@ -50,7 +66,6 @@ double toc() {
 #endif
 }
 
-
 // This is a typedef for a random number generator.
 // Try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand
 typedef boost::minstd_rand _rnd_gen_type;
index 51d026e..49b5352 100644 (file)
@@ -24,7 +24,7 @@ CHANGES:
 
 #include <cmath>
 #include <cstring>
-
+#include <dai/util.h>
 
 
 namespace x2x {
index 02ca982..d33880f 100644 (file)
@@ -25,6 +25,7 @@
 #include <numeric>
 #include <cmath>
 #include <cstdlib>
+#include <cstring>
 #include <boost/program_options.hpp>
 #include <dai/util.h>
 #include <dai/alldai.h>
@@ -46,8 +47,9 @@ class TestAI {
         double          logZ;
         double          maxdiff;
         double          time;
+        bool            has_logZ;
 
-        TestAI( const FactorGraph &fg, const string &_name, const PropertySet &opts ) : obj(NULL), name(_name), err(), q(), logZ(0.0), maxdiff(0.0), time(0) {
+        TestAI( const FactorGraph &fg, const string &_name, const PropertySet &opts ) : obj(NULL), name(_name), err(), q(), logZ(0.0), maxdiff(0.0), time(0), has_logZ(true) {
             double tic = toc();
             obj = newInfAlg( name, fg, opts );
             time += toc() - tic;
@@ -90,8 +92,8 @@ class TestAI {
 
         vector<Factor> allBeliefs() {
             vector<Factor> result;
-            for( size_t i = 0; i < obj->nrVars(); i++ )
-                result.push_back( obj->belief( obj->var(i) ) );
+            for( size_t i = 0; i < obj->fg().nrVars(); i++ )
+                result.push_back( obj->belief( obj->fg().var(i) ) );
             return result;
         }
 
@@ -104,7 +106,12 @@ class TestAI {
                 obj->init();
                 obj->run();
                 time += toc() - tic;
-                logZ = real(obj->logZ());
+                try {
+                    logZ = obj->logZ();
+                    has_logZ = true;
+                } catch( Exception &e ) {
+                    has_logZ = false;
+                }
                 maxdiff = obj->maxDiff();
                 q = allBeliefs();
             };
@@ -148,10 +155,10 @@ pair<string, PropertySet> parseMethod( const string &_s, const map<string,string
     if( pos == string::npos )
         throw "Malformed method";
     size_t n = 0;
-    for( ; n < sizeof(DAINames) / sizeof(string); n++ )
+    for( ; strlen( DAINames[n] ) != 0; n++ )
         if( name == DAINames[n] )
             break;
-    if( n == sizeof(DAINames) / sizeof(string) )
+    if( strlen( DAINames[n] ) == 0 )
         throw "Unknown inference algorithm";
 
     stringstream ss;
@@ -178,6 +185,7 @@ int main( int argc, char *argv[] ) {
         double tol;
         size_t maxiter;
         size_t verbose;
+        bool report_time = true;
 
         po::options_description opts_required("Required options");
         opts_required.add_options()
@@ -192,6 +200,7 @@ int main( int argc, char *argv[] ) {
             ("tol", po::value< double >(&tol), "Override tolerance")
             ("maxiter", po::value< size_t >(&maxiter), "Override maximum number of iterations")
             ("verbose", po::value< size_t >(&verbose), "Override verbosity")
+            ("report-time", po::value< bool >(&report_time), "Report calculation time")
         ;
 
         po::options_description cmdline_options;
@@ -203,7 +212,7 @@ int main( int argc, char *argv[] ) {
 
         if( vm.count("help") || !(vm.count("filename") && vm.count("methods")) ) {
             cout << "Reads factorgraph <filename.fg> and performs the approximate" << endl;
-            cout << "inference algorithms <method*>, reporting clocks, max and average" << endl;
+            cout << "inference algorithms <method*>, reporting calculation time, max and average" << endl;
             cout << "error and relative logZ error (comparing with the results of" << endl;
             cout << "<method0>, the base method, for which one can use JTREE_HUGIN)." << endl << endl;
             cout << opts_required << opts_optional << endl;
@@ -250,8 +259,10 @@ int main( int argc, char *argv[] ) {
             cout << "# " << filename << endl;
             cout.width( 40 );
             cout << left << "# METHOD" << "  ";
-            cout.width( 10 );
-            cout << right << "SECONDS" << "   ";
+            if( report_time ) {
+                cout.width( 10 );
+                cout << right << "SECONDS" << "   ";
+            }
             cout.width( 10 );
             cout << "MAX ERROR" << "  ";
             cout.width( 10 );
@@ -281,8 +292,10 @@ int main( int argc, char *argv[] ) {
                 cout.width( 40 );
 //                cout << left << piet.identify() << "  ";
                 cout << left << methods[m] << "  ";
-                cout.width( 10 );
-                cout << right << piet.time << "    ";
+                if( report_time ) {
+                    cout.width( 10 );
+                    cout << right << piet.time << "    ";
+                }
 
                 if( m > 0 ) {
                     cout.setf( ios_base::scientific );
@@ -294,8 +307,12 @@ int main( int argc, char *argv[] ) {
                     double ae = clipdouble( piet.avgErr(), 1e-9 );
                     cout << ae << "  ";
                     cout.width( 10 );
-                    double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
-                    cout << le << "  ";
+                    if( piet.has_logZ ) {
+                        double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
+                        cout << le << "  ";
+                    } else {
+                        cout << "N/A         ";
+                    }
                     cout.width( 10 );
                     double md = clipdouble( piet.maxdiff, 1e-9 );
                     if( isnan( me ) )
index 2cf6099..aecb684 100755 (executable)
@@ -1,2 +1,2 @@
 #!/bin/bash
-./test --aliases aliases.conf --filename $1 --methods JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND
+./test --report-time false --aliases aliases.conf --filename $1 --methods JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND
diff --git a/tests/testall.bat b/tests/testall.bat
new file mode 100755 (executable)
index 0000000..dc568a9
--- /dev/null
@@ -0,0 +1 @@
+test --report-time false --aliases aliases.conf --filename %1 --methods JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND
index ac36f23..ff7e880 100644 (file)
@@ -1,43 +1,43 @@
 # testfast.fg
-# METHOD                                      CLOCKS     MAX ERROR   AVG ERROR  LOGZ ERROR     MAXDIFF
-JTREE_HUGIN                                        2    
-JTREE_SHSH                                         2     1.000e-09   1.000e-09   1.000e-09   1.000e-09
-BP_SEQFIX                                          1     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_SEQRND                                          1     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_SEQMAX                                          0     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_PARALL                                          1     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_SEQFIX_LOG                                      1     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_SEQRND_LOG                                      0     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_SEQMAX_LOG                                      1     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-BP_PARALL_LOG                                      1     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-MF_SEQRND                                          3     3.607e-01   1.904e-01  -9.409e-02   1.000e-09
-TREEEP                                             5     3.268e-02   8.023e-03   6.340e-04   1.000e-09
-TREEEPWC                                           4     2.356e-02   1.026e-02   4.876e-03   1.000e-09
-GBP_MIN                                            2     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-GBP_DELTA                                          8     6.291e-01   3.350e-01  -3.303e-01   1.000e-09
-GBP_LOOP3                                          2     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-GBP_LOOP4                                          6     7.893e-01   3.569e-01  -3.781e-01   1.000e-09
-HAK_MIN                                           27     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-HAK_DELTA                                        394     3.684e-01   1.892e-01   9.675e-01   1.000e-09
-HAK_LOOP3                                         27     9.483e-02   3.078e-02   1.737e-02   1.000e-09
-HAK_LOOP4                                        726     4.970e-03   1.486e-03  -2.503e-04   1.000e-09
-MR_RESPPROP_FULL                                  36     1.676e-02   4.933e-03         nan   1.000e-09
-MR_CLAMPING_FULL                                  22     8.359e-03   2.508e-03         nan   1.000e-09
-MR_EXACT_FULL                                     20     3.527e-03   1.038e-03         nan   1.000e-09
-MR_RESPPROP_LINEAR                                37     1.932e-02   5.506e-03         nan   1.000e-09
-MR_CLAMPING_LINEAR                                21     1.089e-02   3.076e-03         nan   1.000e-09
-MR_EXACT_LINEAR                                   20     5.617e-03   1.742e-03         nan   1.000e-09
-LCBP_FULLCAV_SEQFIX                               31     1.225e-03   5.589e-04         nan   1.000e-09
-LCBP_FULLCAVin_SEQFIX                             33     1.225e-03   5.589e-04         nan   1.000e-09
-LCBP_FULLCAV_SEQRND                               30     1.225e-03   5.589e-04         nan   1.000e-09
-LCBP_FULLCAVin_SEQRND                             33     1.225e-03   5.589e-04         nan   1.000e-09
-LCBP_FULLCAV_NONE                                 25     1.318e-02   2.644e-03         nan   1.000e+00
-LCBP_FULLCAVin_NONE                               26     1.318e-02   2.644e-03         nan   1.000e+00
-LCBP_PAIRCAV_SEQFIX                               27     1.564e-02   5.284e-03         nan   1.000e-09
-LCBP_PAIRCAVin_SEQFIX                             27     1.564e-02   5.284e-03         nan   1.000e-09
-LCBP_PAIRCAV_SEQRND                               25     1.564e-02   5.284e-03         nan   1.000e-09
-LCBP_PAIRCAVin_SEQRND                             26     1.564e-02   5.284e-03         nan   1.000e-09
-LCBP_PAIRCAV_NONE                                 19     1.869e-01   6.816e-02         nan   1.000e+00
-LCBP_PAIRCAVin_NONE                               20     1.869e-01   6.816e-02         nan   1.000e+00
-LCBP_UNICAV_SEQFIX                                 6     9.483e-02   3.078e-02         nan   1.000e-09
-LCBP_UNICAV_SEQRND                                 8     9.483e-02   3.078e-02         nan   1.000e-09
+# METHOD                                  MAX ERROR   AVG ERROR   LOGZ ERROR  MAXDIFF   
+JTREE_HUGIN                               
+JTREE_SHSH                                1.000e-09   1.000e-09   1.000e-09   1.000e-09 
+BP_SEQFIX                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_SEQRND                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_SEQMAX                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_PARALL                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_SEQFIX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_SEQRND_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_SEQMAX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+BP_PARALL_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+MF_SEQRND                                 3.607e-01   1.904e-01   -9.409e-02  1.000e-09 
+TREEEP                                    3.268e-02   8.023e-03   6.340e-04   1.000e-09 
+TREEEPWC                                  2.356e-02   1.026e-02   4.876e-03   1.000e-09 
+GBP_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+GBP_DELTA                                 6.291e-01   3.350e-01   -3.303e-01  1.000e-09 
+GBP_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+GBP_LOOP4                                 7.893e-01   3.569e-01   -3.781e-01  1.000e-09 
+HAK_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+HAK_DELTA                                 3.684e-01   1.892e-01   9.675e-01   1.000e-09 
+HAK_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09 
+HAK_LOOP4                                 4.970e-03   1.486e-03   -2.503e-04  1.000e-09 
+MR_RESPPROP_FULL                          1.676e-02   4.933e-03   N/A         1.000e-09 
+MR_CLAMPING_FULL                          8.359e-03   2.508e-03   N/A         1.000e-09 
+MR_EXACT_FULL                             3.527e-03   1.038e-03   N/A         1.000e-09 
+MR_RESPPROP_LINEAR                        1.932e-02   5.506e-03   N/A         1.000e-09 
+MR_CLAMPING_LINEAR                        1.089e-02   3.076e-03   N/A         1.000e-09 
+MR_EXACT_LINEAR                           5.617e-03   1.742e-03   N/A         1.000e-09 
+LCBP_FULLCAV_SEQFIX                       1.225e-03   5.589e-04   N/A         1.000e-09 
+LCBP_FULLCAVin_SEQFIX                     1.225e-03   5.589e-04   N/A         1.000e-09 
+LCBP_FULLCAV_SEQRND                       1.225e-03   5.589e-04   N/A         1.000e-09 
+LCBP_FULLCAVin_SEQRND                     1.225e-03   5.589e-04   N/A         1.000e-09 
+LCBP_FULLCAV_NONE                         1.318e-02   2.644e-03   N/A         1.000e+00 
+LCBP_FULLCAVin_NONE                       1.318e-02   2.644e-03   N/A         1.000e+00 
+LCBP_PAIRCAV_SEQFIX                       1.564e-02   5.284e-03   N/A         1.000e-09 
+LCBP_PAIRCAVin_SEQFIX                     1.564e-02   5.284e-03   N/A         1.000e-09 
+LCBP_PAIRCAV_SEQRND                       1.564e-02   5.284e-03   N/A         1.000e-09 
+LCBP_PAIRCAVin_SEQRND                     1.564e-02   5.284e-03   N/A         1.000e-09 
+LCBP_PAIRCAV_NONE                         1.869e-01   6.816e-02   N/A         1.000e+00 
+LCBP_PAIRCAVin_NONE                       1.869e-01   6.816e-02   N/A         1.000e+00 
+LCBP_UNICAV_SEQFIX                        9.483e-02   3.078e-02   N/A         1.000e-09 
+LCBP_UNICAV_SEQRND                        9.483e-02   3.078e-02   N/A         1.000e-09 
index 22f71bb..55630a1 100755 (executable)
@@ -1,16 +1,8 @@
 #!/bin/bash
 TMPFILE1=`mktemp /var/tmp/testfast.XXXXXX`
 trap 'rm -f $TMP_FILE' 0 1 15
-TMPFILE2=`mktemp /var/tmp/testfast.XXXXXX`
-trap 'rm -f $TMP_FILE2' 0 1 15
-TMPFILE3=`mktemp /var/tmp/testfast.XXXXXX`
-trap 'rm -f $TMP_FILE3' 0 1 15
 
 ./testall testfast.fg > $TMPFILE1
-awk '/^#/ {} NF == 2 && !/^#/ {} NF == 6 && !/^#/ { print $1, $3, $4, $5, $6; }' $TMPFILE1 > $TMPFILE2
-awk '/^#/ {} NF == 2 && !/^#/ {} NF == 6 && !/^#/ { print $1, $3, $4, $5, $6; }' testfast.out > $TMPFILE3
-diff -s $TMPFILE2 $TMPFILE3
+diff -s $TMPFILE1 testfast.out
 
 rm -f $TMPFILE1
-rm -f $TMPFILE2
-rm -f $TMPFILE3
index cc5fe80..4feb946 100644 (file)
 
 
 #include <iostream>
+#include <fstream>
+#include <vector>
 #include <iterator>
+#include <algorithm>
 #include <boost/program_options.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/io.hpp>
 #include <dai/factorgraph.h>
 #include <dai/weightedgraph.h>
 #include <dai/util.h>
 using namespace std;
 using namespace dai;
 namespace po = boost::program_options;
+typedef boost::numeric::ublas::compressed_matrix<double> matrix;
+typedef matrix::value_array_type::const_iterator matrix_vcit;
+typedef matrix::index_array_type::const_iterator matrix_icit;
+
+
+Factor BinaryFactor( const Var &n, double field ) {
+    assert( n.states() == 2 );
+    double buf[2];
+    buf[0] = exp(-field); 
+    buf[1] = exp(field);
+    return Factor(n, &buf[0]);
+}
+
+
+Factor BinaryFactor( const Var &n1, const Var &n2, double coupling ) {
+    assert( n1.states() == 2 );
+    assert( n2.states() == 2 );
+    assert( n1 != n2 );
+    double buf[4];
+    buf[0] = (buf[3] = exp(coupling));
+    buf[1] = (buf[2] = exp(-coupling));
+    return Factor(n1 | n2, &buf[0]);
+}
+
+
+Factor RandomFactor( const VarSet &ns, double beta ) {
+    Factor fac( ns );
+    for( size_t t = 0; t < fac.states(); t++ )
+        fac[t] = exp(rnd_stdnormal() * beta);
+    return fac;
+}
+
+
+Factor PottsFactor( const Var &n1, const Var &n2, double beta ) {
+    Factor fac( n1 | n2, 1.0 );
+    assert( n1.states() == n2.states() );
+    for( size_t s = 0; s < n1.states(); s++ )
+        fac[ s * (n1.states() + 1) ] = exp(beta);
+    return fac;
+}
 
 
 void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) {
        vector<Var> vars;
        vector<Factor> factors;
 
+    vars.reserve(N);
        for( size_t i = 0; i < N; i++ )
                vars.push_back(Var(i,2));
 
@@ -45,286 +91,689 @@ void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) {
                        do {
                                size_t newind = (size_t)(N * rnd_uniform());
                                Var newvar = Var(newind, 2);
-                               if( !(vars.contains( newvar )) ) {
+                               if( !vars.contains( newvar ) ) {
                                        vars |= newvar;
                                        break;
                                }
                        } while( 1 );
                }
-               Factor newfac(vars);
-               for( size_t t = 0; t < newfac.states(); t++ )
-                       newfac[t] = exp(rnd_stdnormal() * sigma);
-               factors.push_back(newfac);
+        factors.push_back( RandomFactor( vars, sigma ) );
        }
 
-    fg = FactorGraph(factors);
-};
+    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+}
 
 
-void MakeFullFG( size_t N, double sigma_w, double sigma_th, string type, FactorGraph &fg ) {
-    vector<Var> vars;
+// w should be upper triangular or lower triangular
+void WTh2FG( const matrix &w, const vector<double> &th, FactorGraph &fg ) {
+    vector<Var>    vars;
     vector<Factor> factors;
 
-    double w[N][N];
-    double th[N];
-    double buf[4];
+    size_t N = th.size();
+    assert( (w.size1() == N) && (w.size2() == N) );
 
+    vars.reserve(N);
     for( size_t i = 0; i < N; i++ )
         vars.push_back(Var(i,2));
-    
-    for( size_t i = 0; i < N; i++ )
-        for( size_t j = 0; j < N; j++ )
-            w[i][j] = 0.0;
 
+    factors.reserve( w.nnz() + N );
+    // walk through the sparse array structure
+    // this is similar to matlab sparse arrays
+    // index2 gives the column index (similar to ir in matlab)
+    // index1 gives the starting indices for each row (similar to jc in matlab)
+    size_t i = 0;
+    for( size_t pos = 0; pos < w.nnz(); pos++ ) {
+        while( pos == w.index1_data()[i+1] )
+            i++;
+        size_t j = w.index2_data()[pos];
+        double w_ij = w.value_data()[pos];
+        factors.push_back( BinaryFactor( vars[i], vars[j], w_ij ) );
+    }
     for( size_t i = 0; i < N; i++ )
-        for( size_t j = i+1; j < N; j++ ) {
-            w[i][j] = rnd_stdnormal() * sigma_w;
-            if( type == "fe" )
-                w[i][j] = fabs(w[i][j]);
-            else if( type == "af" )
-                w[i][j] = -fabs(w[i][j]);
-            w[j][i] = w[i][j];
-            buf[0] = (buf[3] = exp(w[i][j]));
-            buf[1] = (buf[2] = exp(-w[i][j]));
-            factors.push_back(Factor(VarSet(vars[i],vars[j]),buf));
-        }
+        factors.push_back( BinaryFactor( vars[i], th[i] ) );
+
+    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+}
+
+
+void MakeFullFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
+    matrix w(N,N,N*(N-1)/2);
+    vector<double> th(N,0.0);
 
     for( size_t i = 0; i < N; i++ ) {
-        th[i] = rnd_stdnormal() * sigma_th;
-        buf[0] = exp(th[i]);
-        buf[1] = exp(-th[i]);
-        factors.push_back(Factor(vars[i],buf));
+        for( size_t j = i+1; j < N; j++ )
+            w(i,j) = rnd_stdnormal() * sigma_w + mean_w;
+        th[i] = rnd_stdnormal() * sigma_th + mean_th;
     }
 
-    fg = FactorGraph(factors);
-};
+    WTh2FG( w, th, fg );
+}
 
 
-void MakeGridFG( long periodic, long n, double sigma_w, double sigma_th, string type, FactorGraph &fg ) {
+void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, double beta, FactorGraph &fg ) {
     vector<Var> vars;
     vector<Factor> factors;
+    
+    for( size_t i1 = 0; i1 < n1; i1++ )
+        for( size_t i2 = 0; i2 < n2; i2++ )
+            for( size_t i3 = 0; i3 < n3; i3++ ) {
+                vars.push_back( Var( i1*n2*n3 + i2*n3 + i3, states ) );
+                if( i1 )
+                    factors.push_back( PottsFactor( vars.back(), vars[ (i1-1)*n2*n3 + i2*n3 + i3 ], beta ) );
+                if( i2 )
+                    factors.push_back( PottsFactor( vars.back(), vars[ i1*n2*n3 + (i2-1)*n3 + i3 ], beta ) );
+                if( i3 )
+                    factors.push_back( PottsFactor( vars.back(), vars[ i1*n2*n3 + i2*n3 + (i3-1) ], beta ) );
+            }
 
-    long N = n*n;
+    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+}
 
-    double w[N][N];
-    double th[N];
-    double buf[4];
 
-    for( long i = 0; i < N; i++ )
-        vars.push_back(Var(i,2));
-    
-    for( long i = 0; i < N; i++ )
-        for( long j = 0; j < N; j++ )
-            w[i][j] = 0.0;
+void MakeGridFG( long periodic, size_t n, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
+    size_t N = n*n;
 
-    for( long i = 0; i < n; i++ )
-        for( long j = 0; j < n; j++ ) {
+    matrix w(N,N,2*N);
+    vector<double> th(N,0.0);
+
+    for( size_t i = 0; i < n; i++ )
+        for( size_t j = 0; j < n; j++ ) {
             if( i+1 < n || periodic )
-                w[i*n+j][((i+1)%n)*n+j] = 1.0;
-            if( i > 0 || periodic )
-                w[i*n+j][((i+n-1)%n)*n+j] = 1.0;
+                w(i*n+j, ((i+1)%n)*n+j) = rnd_stdnormal() * sigma_w + mean_w;
             if( j+1 < n || periodic )
-                w[i*n+j][i*n+((j+1)%n)] = 1.0;
-            if( j > 0 || periodic )
-                w[i*n+j][i*n+((j+n-1)%n)] = 1.0;
+                w(i*n+j, i*n+((j+1)%n)) = rnd_stdnormal() * sigma_w + mean_w;
+            th[i*n+j] = rnd_stdnormal() * sigma_th + mean_th;
         }
 
-    for( long i = 0; i < N; i++ )
-        for( long j = i+1; j < N; j++ ) 
-            if( w[i][j] ) {
-                w[i][j] = rnd_stdnormal() * sigma_w;
-                if( type == "fe" )
-                    w[i][j] = fabs(w[i][j]);
-                else if( type == "af" )
-                    w[i][j] = -fabs(w[i][j]);
-                w[j][i] = w[i][j];
-                buf[0] = (buf[3] = exp(w[i][j]));
-                buf[1] = (buf[2] = exp(-w[i][j]));
-                factors.push_back(Factor(VarSet(vars[i],vars[j]),buf));
-            }
+    WTh2FG( w, th, fg );
+}
+
+            
+void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, double beta, FactorGraph &fg ) {
+    size_t N = n*n;
+
+    vector<Var>    vars;
+    vector<Factor> factors;
+
+    vars.reserve(N);
+    for( size_t i = 0; i < N; i++ )
+        vars.push_back(Var(i, states));
+
+    factors.reserve( 2 * N );
+    for( size_t i = 0; i < n; i++ ) {
+        for( size_t j = 0; j < n; j++ ) {
+            if( i+1 < n || periodic )
+                factors.push_back( RandomFactor( vars[i*n+j] | vars[((i+1)%n)*n+j], beta ) );
+            if( j+1 < n || periodic )
+                factors.push_back( RandomFactor( vars[i*n+j] | vars[i*n+((j+1)%n)], beta ) );
+        }
+    }
+
+    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+}
 
-    for( long i = 0; i < N; i++ ) {
-        th[i] = rnd_stdnormal() * sigma_th;
-        buf[0] = exp(th[i]);
-        buf[1] = exp(-th[i]);
-        factors.push_back(Factor(vars[i],buf));
+
+void MakeLoopFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
+    matrix w(N,N,N);
+    vector<double> th(N,0.0);
+
+    for( size_t i = 0; i < N; i++ ) {
+        w(i, (i+1)%N) = rnd_stdnormal() * sigma_w + mean_w;
+        th[i] = rnd_stdnormal() * sigma_th + mean_th;
     }
 
-    fg = FactorGraph(factors);
-};
+    WTh2FG( w, th, fg );
+}
 
 
-void MakeDRegFG( size_t N, size_t d, double sigma_w, double sigma_th, string type, FactorGraph &fg ) {
-    vector<Var> vars;
+void MakeLoopNonbinaryFG( size_t N, size_t states, double beta, FactorGraph &fg ) {
+    vector<Var>    vars;
     vector<Factor> factors;
 
-    double w[N][N];
-    double th[N];
-    double buf[4];
-
-    for( size_t i = 0; i < N; i++ )
-        vars.push_back(Var(i,2));
-    
+    vars.reserve(N);
     for( size_t i = 0; i < N; i++ )
-        for( size_t j = 0; j < N; j++ )
-            w[i][j] = 0.0;
+        vars.push_back(Var(i, states));
 
-    UEdgeVec g = RandomDRegularGraph( N, d );
-    for( size_t i = 0; i < g.size(); i++ ) {
-        w[g[i].n1][g[i].n2] = 1.0;
-        w[g[i].n2][g[i].n1] = 1.0;
+    factors.reserve( N );
+    for( size_t i = 0; i < N; i++ ) {
+        factors.push_back( RandomFactor( vars[i] | vars[(i+1)%N], beta ) );
     }
-    
-    for( size_t i = 0; i < N; i++ )
-        for( size_t j = i+1; j < N; j++ ) 
-            if( w[i][j] ) {
-                w[i][j] = rnd_stdnormal() * sigma_w;
-                if( type == "fe" )
-                    w[i][j] = fabs(w[i][j]);
-                else if( type == "af" )
-                    w[i][j] = -fabs(w[i][j]);
-                w[j][i] = w[i][j];
-                buf[0] = (buf[3] = exp(w[i][j]));
-                buf[1] = (buf[2] = exp(-w[i][j]));
-                factors.push_back(Factor(VarSet(vars[i],vars[j]),buf));
-            }
+
+    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+}
+
+
+void MakeTreeFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
+    matrix w(N,N,N-1);
+    vector<double> th(N,0.0);
 
     for( size_t i = 0; i < N; i++ ) {
-        th[i] = rnd_stdnormal() * sigma_th;
-        buf[0] = exp(th[i]);
-        buf[1] = exp(-th[i]);
-        factors.push_back(Factor(vars[i],buf));
+        th[i] = rnd_stdnormal() * sigma_th + mean_th;
+        if( i > 0 ) {
+            size_t j = rnd_int( 0, i-1 );
+            w(i,j) = rnd_stdnormal() * sigma_w + mean_w;
+        }
     }
 
-    fg = FactorGraph(factors);
-};
+    WTh2FG( w, th, fg );
+}
+
+
+void MakeDRegFG( size_t N, size_t d, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
+    matrix w(N,N,(d*N)/2);
+    vector<double> th(N,0.0);
+
+    UEdgeVec g = RandomDRegularGraph( N, d );
+    for( size_t i = 0; i < g.size(); i++ )
+        w(g[i].n1, g[i].n2) = rnd_stdnormal() * sigma_w + mean_w;
+    
+    for( size_t i = 0; i < N; i++ )
+        th[i] = rnd_stdnormal() * sigma_th + mean_th;
+
+    WTh2FG( w, th, fg );
+}
+
 
+// N = number of variables
+// n = size of variable neighborhoods
+// K = number of factors
+// k = size of factor neighborhoods
+// asserts: N * n == K * k
+BipartiteGraph CreateRandomBipartiteGraph( size_t N, size_t K, size_t n, size_t k ) {
+    BipartiteGraph G;
 
-const char *HOITYPE = "hoi";
-const char *FULLTYPE = "full";
-const char *GRIDTYPE = "grid";
-const char *DREGTYPE = "dreg";
+    assert( N * n == K * k );
 
+    // build lists of degree-repeated vertex numbers
+    std::vector<size_t> stubs1(N*n,0);
+    for( size_t i = 0; i < N; i++ )
+        for( size_t t = 0; t < n; t++ )
+            stubs1[i*n + t] = i;
+
+    // build lists of degree-repeated vertex numbers
+    std::vector<size_t> stubs2(K*k,0);
+    for( size_t I = 0; I < K; I++ )
+        for( size_t t = 0; t < k; t++ )
+            stubs2[I*k + t] = I;
+
+    // shuffle lists
+    random_shuffle( stubs1.begin(), stubs1.end() );
+    random_shuffle( stubs2.begin(), stubs2.end() );
 
-// Old usages:
-// create_full_fg <N> <sigma_w> <sigma_th> <subtype>
-// create_grid_fg <periodic> <n> <sigma_w> <sigma_th> <subtype>
-// create_dreg_fg <d> <N> <sigma_w> <sigma_th> <subtype>
+    // add edges
+    vector<BipartiteGraph::Edge> edges;
+    edges.reserve( N*n );
+    for( size_t e = 0; e < N*n; e++ )
+        edges.push_back( BipartiteGraph::Edge(stubs1[e], stubs2[e]) );
+
+    // finish construction
+    G.create( N, K, edges.begin(), edges.end() );
+
+    return G;
+}
+
+
+// Returns x**n % p, assuming p is prime
+int powmod (int x, int n, int p) {
+    int y = 1;
+    for( int m = 0; m < n; m++ )
+        y = (x * y) % p;
+    return y;
+}
+
+
+// Returns order of x in GF(p) with p prime
+size_t order (int x, int p) {
+    x = x % p;
+    assert( x != 0 );
+    size_t n = 0;
+    size_t prod = 1;
+    do {
+        prod = (prod * x) % p;
+        n++;
+    } while( prod != 1 ); 
+    return n;
+}
+
+
+// Returns whether n is a prime number
+bool isPrime (size_t n) {
+    bool result = true;
+    for( size_t k = 2; (k < n) && result; k++ )
+        if( n % k == 0 )
+            result = false;
+    return result;
+}
+
+
+// Make a regular LDPC graph with N=6, j=2, K=4, k=3
+BipartiteGraph CreateSmallLDPCGraph() {
+    BipartiteGraph G;
+    size_t N=4, j=3, K=4; // k=3;
+
+    typedef BipartiteGraph::Edge Edge;
+    vector<Edge> edges;
+    edges.reserve( N*j );
+    edges.push_back( Edge(0,0) ); edges.push_back( Edge(1,0) ); edges.push_back( Edge(2,0) );
+    edges.push_back( Edge(0,1) ); edges.push_back( Edge(1,1) ); edges.push_back( Edge(3,1) );
+    edges.push_back( Edge(0,2) ); edges.push_back( Edge(2,2) ); edges.push_back( Edge(3,2) );
+    edges.push_back( Edge(1,3) ); edges.push_back( Edge(2,3) ); edges.push_back( Edge(3,3) );
+
+    // finish construction
+    G.create( N, K, edges.begin(), edges.end() );
+
+    return G;
+}
+
+
+// Use construction described in "A Class of Group-Structured LDPC Codes"
+// by R. M. Tanner, D. Sridhara and T. Fuja
+// Proceedings of ICSTA, 2001
+//
+// Example parameters: (p,j,k) = (31,3,5)
+// j and k must be divisors of p-1
+BipartiteGraph CreateGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) {
+    BipartiteGraph G;
+
+    size_t n = j;
+    size_t N = p * k;
+    size_t K = p * j;
+
+    size_t a, b;
+    for( a = 2; a < p; a++ )
+        if( order(a,p) == k )
+            break;
+    assert( a != p );
+    for( b = 2; b < p; b++ )
+        if( order(b,p) == j )
+            break;
+    assert( b != p );
+//    cout << "# order(a=" << a << ") = " << order(a,p) << endl;
+//    cout << "# order(b=" << b << ") = " << order(b,p) << endl;
+
+    assert( N * n == K * k );
+
+    typedef BipartiteGraph::Edge Edge;
+    vector<Edge> edges;
+    edges.reserve( N * n );
+
+    for( size_t s = 0; s < j; s++ )
+        for( size_t t = 0; t < k; t++ ) {
+            size_t P = (powmod(b,s,p) * powmod(a,t,p)) % p;
+            for( size_t m = 0; m < p; m++ )
+                edges.push_back( Edge(t*p + m, s*p + ((m + P) % p)) );
+        }
+
+    // finish construction
+    G.create( N, K, edges.begin(), edges.end() );
+
+    return G;
+}
+
+
+// Make parity check table
+void MakeParityCheck( double *result, size_t n, double eps ) {
+    size_t N = 1 << n;
+    for( size_t i = 0; i < N; i++ ) {
+        size_t c = 0;
+        for( size_t t = 0; t < n; t++ )
+            if( i & (1 << t) )
+                c ^= 1;
+        if( c )
+            result[i] = eps;
+        else
+            result[i] = 1.0 - eps;
+    }
+    return;
+}
+
+
+const char *FULL_TYPE = "full";
+const char *GRID_TYPE = "grid";
+const char *GRID_TORUS_TYPE = "grid_torus";
+const char *DREG_TYPE = "dreg";
+const char *HOI_TYPE = "hoi";
+const char *LDPC_RANDOM_TYPE = "ldpc_random";
+const char *LDPC_GROUP_TYPE = "ldpc_group";
+const char *LDPC_SMALL_TYPE = "ldpc_small";
+const char *LOOP_TYPE = "loop";
+const char *TREE_TYPE = "tree";
+const char *POTTS3D_TYPE = "potts3d";
 
 
 int main( int argc, char *argv[] ) {
     try {
-               size_t N, M, k, d;
-        size_t periodic;
+               size_t N, K, k, d, j, n1, n2, n3;
+        size_t prime;
         size_t seed;
-               double beta, sigma_w, sigma_th;
-        string type, subtype;
+               double beta, sigma_w, sigma_th, noise, mean_w, mean_th;
+        string type;
+        size_t states = 2;
 
         // Declare the supported options.
         po::options_description desc("Allowed options");
         desc.add_options()
             ("help",     "produce help message")
-            ("type",     po::value<string>(&type),     "factor graph type:\n\t'full', 'grid', 'dreg' or 'hoi'")
-            ("seed",     po::value<size_t>(&seed),     "random number seed")
-            ("subtype",  po::value<string>(&subtype),  "interactions type:\n\t'sg', 'fe' or 'af'\n\t(ignored for type=='hoi')")
-            ("N",        po::value<size_t>(&N),        "number of (binary) variables")
-            ("M",        po::value<size_t>(&M),        "number of factors\n\t(only for type=='hoi')")
-            ("k",        po::value<size_t>(&k),        "connectivity of the factors\n\t(only for type=='hoi')")
+            ("type",     po::value<string>(&type),     "factor graph type:\n\t'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d'")
+            ("seed",     po::value<size_t>(&seed),     "random number seed (tries to read from /dev/urandom if not specified)")
+            ("N",        po::value<size_t>(&N),        "number of variables (not for type=='ldpc_small')")
+            ("n1",       po::value<size_t>(&n1),       "width of 3D grid (only for type=='potts3d')")
+            ("n2",       po::value<size_t>(&n2),       "height of 3D grid (only for type=='potts3d')")
+            ("n3",       po::value<size_t>(&n3),       "length of 3D grid (only for type=='potts3d')")
+            ("K",        po::value<size_t>(&K),        "number of factors\n\t(only for type=='hoi' and 'type=='ldpc_{random,group}')")
+            ("k",        po::value<size_t>(&k),        "number of variables per factor\n\t(only for type=='hoi' and type=='ldpc_{random,group}')")
             ("d",        po::value<size_t>(&d),        "variable connectivity\n\t(only for type=='dreg')")
-            ("beta",     po::value<double>(&beta),     "stddev of log-factor entries\n\t(only for type=='hoi')")
-            ("sigma_w",  po::value<double>(&sigma_w),  "stddev of pairwise interactions w_{ij}\n\t(ignored for type=='hoi')")
-            ("sigma_th", po::value<double>(&sigma_th), "stddev of singleton interactions th_i\n\t(ignored for type=='hoi')")
-            ("periodic", po::value<size_t>(&periodic), "0/1 corresponding to nonperiodic/periodic grid\n\t(only for type=='grid')")
+            ("j",        po::value<size_t>(&j),        "number of parity checks per bit\n\t(only for type=='ldpc_{random,group}')")
+            ("prime",    po::value<size_t>(&prime),    "prime number for construction of LDPC code\n\t(only for type=='ldpc_group')")
+            ("beta",     po::value<double>(&beta),     "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)")
+            ("mean_w",   po::value<double>(&mean_w),   "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
+            ("mean_th",  po::value<double>(&mean_th),  "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
+            ("sigma_w",  po::value<double>(&sigma_w),  "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
+            ("sigma_th", po::value<double>(&sigma_th), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'")
+            ("noise",    po::value<double>(&noise),    "bitflip probability for binary symmetric channel (only for type=='ldpc')")
+            ("states",   po::value<size_t>(&states),   "number of states of each variable (should be 2 for all but type=='grid', 'grid_torus', 'loop', 'potts3d')")
         ;
 
-        po::variables_map vm;        
+        po::variables_map vm;
         po::store(po::parse_command_line(argc, argv, desc), vm);
-        po::notify(vm);    
+        po::notify(vm);
 
         if( vm.count("help") || !vm.count("type") ) {
             if( vm.count("type") ) {
-                if( type == HOITYPE ) {
+                if( type == FULL_TYPE ) {
+                    cout << "Creates fully connected pairwise graphical model of <N> binary variables;" << endl;
+                } else if( type == GRID_TYPE ) {
+                    cout << "Creates (non-periodic) 2D Ising grid of (approx.) <N> variables (which need not be binary);" << endl;
+                } else if( type == GRID_TORUS_TYPE ) {
+                    cout << "Creates periodic 2D Ising grid of (approx.) <N> variables (which need not be binary);" << endl;
+                } else if( type == DREG_TYPE ) {
+                    cout << "Creates random d-regular graph of <N> binary variables with uniform degree <d>" << endl;
+                    cout << "(where <d><N> should be even);" << endl;
+                } else if( type == LOOP_TYPE ) {
+                    cout << "Creates a pairwise graphical model consisting of a single loop of" << endl;
+                    cout << "<N> variables (which need not be binary);" << endl;
+                } else if( type == TREE_TYPE ) {
+                    cout << "Creates a pairwise, connected graphical model without cycles (i.e., a tree)" << endl;
+                    cout << "of <N> binary variables;" << endl;
+                } else if( type == HOI_TYPE ) {
                     cout << "Creates a random factor graph of <N> binary variables and" << endl;
-                    cout << "<M> factors, each factor being an interaction of <k> variables." << endl;
+                    cout << "<K> factors, each factor being an interaction of <k> variables." << endl;
                     cout << "The entries of the factors are exponentials of i.i.d. Gaussian" << endl;
                     cout << "variables with mean 0 and standard deviation <beta>." << endl;
-                } else if( type == FULLTYPE ) {
-                    cout << "Creates fully connected pairwise graphical model of <N> variables;" << endl;
-                } else if( type == GRIDTYPE ) {
-                    cout << "Creates 2D Ising grid (periodic if <periodic>!=0) of (approx.) <N> variables;" << endl;
-                } else if( type == DREGTYPE ) {
-                    cout << "Creates random d-regular graph of <N> nodes with uniform degree <d>" << endl;
-                    cout << "(where <d><N> should be even);" << endl;
+                } else if( type == LDPC_RANDOM_TYPE ) {
+                    cout << "Simulates LDPC decoding problem, using a LDPC code of <N> bits and <K> parity" << endl;
+                    cout << "checks, with <k> bits per check and <j> checks per bit, transmitted on a binary" << endl;
+                    cout << "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl;
+                    cout << "codeword has all bits set to zero. The LDPC code is randomly generated." << endl;
+                } else if( type == LDPC_GROUP_TYPE ) {
+                    cout << "Simulates LDPC decoding problem, using a LDPC code of <N> bits and <K> parity" << endl;
+                    cout << "checks, with <k> bits per check and <j> checks per bit, transmitted on a binary" << endl;
+                    cout << "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl;
+                    cout << "codeword has all bits set to zero. The LDPC code is constructed (using group" << endl;
+                    cout << "theory) using a parameter <prime>; <j> and <k> should both be divisors of <prime>-1." << endl;
+                } else if( type == LDPC_SMALL_TYPE ) {
+                    cout << "Simulates LDPC decoding problem, using a LDPC code of 4 bits and 4 parity" << endl;
+                    cout << "checks, with 3 bits per check and 3 checks per bit, transmitted on a binary" << endl;
+                    cout << "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl;
+                    cout << "codeword has all bits set to zero. The LDPC code is fixed." << endl;
+                } else if( type == POTTS3D_TYPE ) {
+                    cout << "Builds 3D Potts model of size <n1>x<n2>x<n3> with nearest-neighbour Potts" << endl;
+                    cout << "interactions with <states> states and inverse temperature <beta>." << endl;
                 } else
-                    cerr << "Unknown type (should be one of 'full', 'grid', 'dreg' or 'hoi')" << endl;
+                    cerr << "Unknown type (should be one of 'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d')" << endl;
                 
-                if( type == FULLTYPE || type == GRIDTYPE || type == DREGTYPE ) {
-                    cout << "singleton interactions are Gaussian with mean 0 and standard" << endl;
-                    cout << "deviation <sigma_th>; pairwise interactions are Gaussian with mean 0" << endl;
-                    cout << "and standard deviation <sigma_w> if <subtype>=='sg', absolute value" << endl;
-                    cout << "is taken if <subtype>=='fe' and a minus sign is added if <subtype>=='af'." << endl;
-                }
+                if( type == FULL_TYPE || type == GRID_TYPE || type == GRID_TORUS_TYPE || type == DREG_TYPE || type == LOOP_TYPE || type == TREE_TYPE ) {
+                    if( type == GRID_TYPE || type == GRID_TORUS_TYPE || type == LOOP_TYPE ) {
+                        cout << "if <states> > 2: factor entries are exponents of Gaussians with mean 0 and standard deviation beta; otherwise," << endl;
+                    }
+                    cout << "singleton interactions are Gaussian with mean <mean_th> and standard" << endl;
+                    cout << "deviation <sigma_th>; pairwise interactions are Gaussian with mean" << endl;
+                    cout << "<mean_w> and standard deviation <sigma_w>." << endl;
+                } 
             }
             cout << endl << desc << endl;
             return 1;
         }
 
-        if( !vm.count("seed") )
-            throw "Please specify random number seed.";
+        if( !vm.count("states") )
+            states = 2;
+
+        if( !vm.count("seed") ) {
+            ifstream infile;
+            bool success;
+            infile.open( "/dev/urandom" );
+            success = infile.is_open();
+            if( success ) {
+                infile.read( (char *)&seed, sizeof(size_t) / sizeof(char) );
+                success = infile.good();
+                infile.close();
+            }
+            if( !success )
+                throw "Please specify random number seed.";
+        }
         rnd_seed( seed );
-//      srand( gsl_rng_default_seed );
 
         FactorGraph fg;
 
         cout << "# Factor graph made by " << argv[0] << endl;
         cout << "# type = " << type << endl;
 
-        if( type == HOITYPE ) {
-            if( !vm.count("N") || !vm.count("M") || !vm.count("k") || !vm.count("beta") )
-                throw "Please specify all required arguments";
-            do {
-                MakeHOIFG( N, M, k, beta, fg );
-            } while( !fg.G.isConnected() );
-
-            cout << "# N = " << N << endl;
-            cout << "# M = " << M << endl;
-            cout << "# k = " << k << endl;
-            cout << "# beta = " << beta << endl;
-        } else if( type == FULLTYPE ) {
-            if( !vm.count("N") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("subtype") )
+        if( type == FULL_TYPE ) {
+            if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
                 throw "Please specify all required arguments";
-            MakeFullFG( N, sigma_w, sigma_th, subtype, fg );
+            MakeFullFG( N, mean_w, mean_th, sigma_w, sigma_th, fg );
 
             cout << "# N = " << N << endl;
+            cout << "# mean_w = " << mean_w << endl;
+            cout << "# mean_th = " << mean_th << endl;
             cout << "# sigma_w = " << sigma_w << endl;
             cout << "# sigma_th = " << sigma_th << endl;
-            cout << "# subtype = " << subtype << endl;
-        } else if( type == GRIDTYPE ) {
-            if( !vm.count("N") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("subtype") || !vm.count("periodic") )
-                throw "Please specify all required arguments";
+        } else if( type == GRID_TYPE || type == GRID_TORUS_TYPE ) {
+            if( states > 2 ) {
+                if( !vm.count("N") || !vm.count("beta") )
+                    throw "Please specify all required arguments";
+            } else {
+                if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
+                    throw "Please specify all required arguments";
+            }
 
             size_t n = (size_t)sqrt((long double)N);
             N = n * n;
 
-            MakeGridFG( periodic, n, sigma_w, sigma_th, subtype, fg );
-            
-            cout << "# periodic = " << periodic << endl;
+            bool periodic = false;
+            if( type == GRID_TYPE )
+                periodic = false;
+            else
+                periodic = true;
+
+            if( states > 2 )
+                MakeGridNonbinaryFG( periodic, n, states, beta, fg );
+            else
+                MakeGridFG( periodic, n, mean_w, mean_th, sigma_w, sigma_th, fg );
+                
             cout << "# n = " << n << endl;
             cout << "# N = " << N << endl;
-            cout << "# sigma_w = " << sigma_w << endl;
-            cout << "# sigma_th = " << sigma_th << endl;
-            cout << "# subtype = " << subtype << endl;
-        } else if( type == DREGTYPE ) {
-            if( !vm.count("N") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("subtype") || !vm.count("d") )
+            
+            if( states > 2 )
+                cout << "# beta = " << beta << endl;
+            else {
+                cout << "# mean_w = " << mean_w << endl;
+                cout << "# mean_th = " << mean_th << endl;
+                cout << "# sigma_w = " << sigma_w << endl;
+                cout << "# sigma_th = " << sigma_th << endl;
+            }
+        } else if( type == DREG_TYPE ) {
+            if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("d") )
                 throw "Please specify all required arguments";
 
-            MakeDRegFG( N, d, sigma_w, sigma_th, subtype, fg );
+            MakeDRegFG( N, d, mean_w, mean_th, sigma_w, sigma_th, fg );
 
             cout << "# N = " << N << endl;
             cout << "# d = " << d << endl;
+            cout << "# mean_w = " << mean_w << endl;
+            cout << "# mean_th = " << mean_th << endl;
+            cout << "# sigma_w = " << sigma_w << endl;
+            cout << "# sigma_th = " << sigma_th << endl;
+        } else if( type == LOOP_TYPE ) {
+            if( states > 2 ) {
+                if( !vm.count("N") || !vm.count("beta") )
+                    throw "Please specify all required arguments";
+            } else {
+                if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
+                    throw "Please specify all required arguments";
+            }
+            if( states > 2 )
+                MakeLoopNonbinaryFG( N, states, beta, fg );
+            else
+                MakeLoopFG( N, mean_w, mean_th, sigma_w, sigma_th, fg );
+
+            cout << "# N = " << N << endl;
+
+            if( states > 2 )
+                cout << "# beta = " << beta << endl;
+            else {
+                cout << "# mean_w = " << mean_w << endl;
+                cout << "# mean_th = " << mean_th << endl;
+                cout << "# sigma_w = " << sigma_w << endl;
+                cout << "# sigma_th = " << sigma_th << endl;
+            }
+        } else if( type == TREE_TYPE ) {
+            if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
+                throw "Please specify all required arguments";
+            MakeTreeFG( N, mean_w, mean_th, sigma_w, sigma_th, fg );
+
+            cout << "# N = " << N << endl;
+            cout << "# mean_w = " << mean_w << endl;
+            cout << "# mean_th = " << mean_th << endl;
             cout << "# sigma_w = " << sigma_w << endl;
             cout << "# sigma_th = " << sigma_th << endl;
-            cout << "# subtype = " << subtype << endl;
+        } else if( type == HOI_TYPE ) {
+            if( !vm.count("N") || !vm.count("K") || !vm.count("k") || !vm.count("beta") )
+                throw "Please specify all required arguments";
+            do {
+                MakeHOIFG( N, K, k, beta, fg );
+            } while( !fg.G.isConnected() );
+
+            cout << "# N = " << N << endl;
+            cout << "# K = " << K << endl;
+            cout << "# k = " << k << endl;
+            cout << "# beta = " << beta << endl;
+        } else if( type == LDPC_RANDOM_TYPE || type == LDPC_GROUP_TYPE || type == LDPC_SMALL_TYPE ) {
+            if( !vm.count("noise") )
+                throw "Please specify all required arguments";
+
+            if( type == LDPC_RANDOM_TYPE ) {
+                if( !vm.count("N") || !vm.count("K") || !vm.count("j") || !vm.count("k") )
+                    throw "Please specify all required arguments";
+
+                if( N * j != K * k )
+                    throw "Parameters should satisfy N * j == K * k";
+            } else if( type == LDPC_GROUP_TYPE ) {
+                if( !vm.count("prime") || !vm.count("j") || !vm.count("k") )
+                    throw "Please specify all required arguments";
+
+                if( !isPrime(prime) )
+                    throw "Parameter <prime> should be prime";
+                if( !((prime-1) % j == 0 ) )
+                    throw "Parameters should satisfy (prime-1) % j == 0";
+                if( !((prime-1) % k == 0 ) )
+                    throw "Parameters should satisfy (prime-1) % k == 0";
+
+                N = prime * k;
+                K = prime * j;
+            } else if( type == LDPC_SMALL_TYPE ) {
+                N = 4;
+                K = 4;
+                j = 3;
+                k = 3;
+            }
+
+            cout << "# N = " << N << endl;
+            cout << "# K = " << K << endl;
+            cout << "# j = " << j << endl;
+            cout << "# k = " << k << endl;
+            if( type == LDPC_GROUP_TYPE )
+                cout << "# prime = " << prime << endl;
+            cout << "# noise = " << noise << endl;
+
+            // p = 31, j = 3, k = 5
+            // p = 37, j = 3, k = 4
+            // p = 7 , j = 2, k = 3
+            // p = 29, j = 2, k = 4
+
+            // Construct likelihood and paritycheck factors
+            double likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise};
+            double *paritycheck = new double[1 << k];
+            MakeParityCheck(paritycheck, k, 0.0);
+
+            // Create LDPC structure
+            BipartiteGraph ldpcG;
+            bool regular;
+            do {
+                if( type == LDPC_GROUP_TYPE )
+                    ldpcG = CreateGroupStructuredLDPCGraph( prime, j, k );
+                else if( type == LDPC_RANDOM_TYPE )
+                    ldpcG = CreateRandomBipartiteGraph( N, K, j, k );
+                else if( type == LDPC_SMALL_TYPE )
+                    ldpcG = CreateSmallLDPCGraph();
+
+                regular = true;
+                for( size_t i = 0; i < N; i++ )
+                    if( ldpcG.nb1(i).size() != j )
+                        regular = false;
+                for( size_t I = 0; I < K; I++ )
+                    if( ldpcG.nb2(I).size() != k )
+                        regular = false;
+            } while( !regular && !ldpcG.isConnected() );
+
+            // Convert to FactorGraph
+            vector<Factor> factors;
+            for( size_t I = 0; I < K; I++ ) {
+                VarSet vs;
+                for( size_t _i = 0; _i < k; _i++ ) {
+                    size_t i = ldpcG.nb2(I)[_i];
+                    vs |= Var( i, 2 );
+                }
+                factors.push_back( Factor( vs, paritycheck ) );
+            }
+            delete paritycheck;
+
+            // Generate noise vector
+            vector<char> noisebits(N,0);
+            size_t bitflips = 0;
+            for( size_t i = 0; i < N; i++ ) {
+                if( rnd_uniform() < noise ) {
+                    noisebits[i] = 1;
+                    bitflips++;
+                }
+            }
+            cout << "# bitflips = " << bitflips << endl;
+
+            // Simulate transmission of all-zero codeword
+            vector<char> input(N,0);
+            vector<char> output(N,0);
+            for( size_t i = 0; i < N; i++ )
+                output[i] = (input[i] + noisebits[i]) & 1;
+
+            // Add likelihoods
+            for( size_t i = 0; i < N; i++ )
+               factors.push_back( Factor(Var(i,2), likelihood + output[i]*2) ); 
+
+            // Construct Factor Graph
+            fg = FactorGraph( factors );
+        } else if( type == POTTS3D_TYPE ) {
+            if( !vm.count("n1") || !vm.count("n2") || !vm.count("n3") || !vm.count("beta") || !vm.count("states") )
+                throw "Please specify all required arguments";
+            Make3DPotts( n1, n2, n3, states, beta, fg );
+
+            cout << "# N = " << n1*n2*n3 << endl;
+            cout << "# n1 = " << n1 << endl;
+            cout << "# n2 = " << n2 << endl;
+            cout << "# n3 = " << n3 << endl;
+            cout << "# beta = " << beta << endl;
+            cout << "# states = " << states << endl;
+        } else {
+            throw "Invalid type";
         }
 
         cout << "# seed = " << seed << endl;