Now compiles also with Visual Studio 2008 under Windows (still buggy!)
authorJoris Mooij <jorism@osun.tuebingen.mpg.de>
Thu, 18 Sep 2008 08:40:52 +0000 (10:40 +0200)
committerJoris Mooij <jorism@osun.tuebingen.mpg.de>
Thu, 18 Sep 2008 08:40:52 +0000 (10:40 +0200)
- Merged utils/createfg from SVN head
- Added Makefile.win for Visual Studio and fixed some bugs;
  HAK doesn't work yet, because of problems with NANs.

Makefile.win [new file with mode: 0755]
include/dai/factorgraph.h
include/dai/lc.h
include/dai/mr.h
include/dai/util.h
include/dai/var.h
src/factorgraph.cpp
src/util.cpp
src/x2x.cpp
tests/test.cpp
utils/createfg.cpp

diff --git a/Makefile.win b/Makefile.win
new file mode 100755 (executable)
index 0000000..7a12f2e
--- /dev/null
@@ -0,0 +1,242 @@
+#   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\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 $(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 $(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 example.exe matlab\*.$(MEXEXT) matlab\*.obj tests/test.exe utils\fg2dot.exe utils\createfg.exe utils\remove_short_loops.exe utils\fginfo.exe $(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
+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\r
+       $(CC) $(CCFLAGS) /Feutils/createfg utils/createfg.cpp factorgraph.obj weightedgraph.obj util.obj bipgraph.obj /link $(LIBS) $(BOOSTLIBS)\r
+\r
+utils/fg2dot.exe : utils/fg2dot.cpp $(HEADERS) factorgraph.obj\r
+       $(CC) $(CCFLAGS) /Feutils/fg2dot utils/fg2dot.cpp factorgraph.obj /link $(LIBS)\r
+\r
+utils/remove_short_loops.exe : utils/remove_short_loops.cpp $(HEADERS) factorgraph.obj\r
+       $(CC) $(CCFLAGS) /Feutils/remove_short_loops utils/remove_short_loops.cpp factorgraph.obj /link $(LIBS)\r
+\r
+utils/fginfo.exe : utils/fginfo.cpp $(HEADERS) factorgraph.obj bipgraph.obj\r
+       $(CC) $(CCFLAGS) /Feutils/fginfo utils/fginfo.cpp factorgraph.obj bipgraph.obj /link $(LIBS)\r
index 77f5f8a..a502ae8 100644 (file)
@@ -25,7 +25,6 @@
 
 #include <iostream>
 #include <map>
-#include <tr1/unordered_map>
 #include <dai/bipgraph.h>
 #include <dai/factor.h>
 
index 231b043..7ded83e 100644 (file)
@@ -81,7 +81,7 @@ 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 { assert( 0 == 1 ); return Factor(); }
         std::vector<Factor> beliefs() const { return _beliefs; }
         Complex logZ() const { return NAN; }
         void CalcBelief (size_t i);
index 153b83a..eb0726d 100644 (file)
@@ -71,7 +71,7 @@ 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 { assert( 0 == 1 ); return Factor(); }
         std::vector<Factor> beliefs() const;
         Complex logZ() const { return NAN; }
         void init() { assert( checkProperties() ); }
@@ -89,7 +89,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); }
 
         bool checkProperties();
 
index d83c204..4c093e0 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 <xmath.h> // for NAN
+    #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 b887316..38a6b5a 100644 (file)
@@ -27,7 +27,6 @@
 #include <string>
 #include <algorithm>
 #include <functional>
-#include <tr1/unordered_map>
 #include <dai/factorgraph.h>
 #include <dai/util.h>
 
@@ -62,7 +61,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;
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 007e7f3..2ac19e4 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>
@@ -148,10 +149,10 @@ pair<string, Properties> 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;
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;