Various changes:
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Wed, 10 Feb 2010 15:45:20 +0000 (16:45 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Wed, 10 Feb 2010 15:45:20 +0000 (16:45 +0100)
* Added GraphAL, an adjacency list implementation for graphs,
  similar to (but simpler than) BipartiteGraph
* Renamed Graph into GraphEL (for Graph, implemented as Edge List)
* Documented and cleaned up utils/createfg
* Renamed some of the factor creating functions introduced in the
  previous commit (now their names start with the prefix 'createFactor')

ChangeLog
Makefile
include/dai/doc.h
include/dai/factor.h
include/dai/graph.h [new file with mode: 0644]
include/dai/weightedgraph.h
src/factor.cpp
src/graph.cpp [new file with mode: 0644]
src/jtree.cpp
src/weightedgraph.cpp
utils/createfg.cpp

index 0c058e6..f18a64a 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,9 @@
-* Added some functionality to create various standard factors to factor.cpp
+* Documented and cleaned up utils/createfg
+* Added GraphAL, an adjacency list implementation for graphs,
+  similar to (but simpler than) BipartiteGraph
+* Renamed Graph into GraphEL (for Graph, implemented as Edge List)
+* Added some functionality to create various standard factors
+  (the functions whose names start with "createFactor")
 * Renamed "Clamped BP" into "Conditioned BP"
 * Improved documentation
 * Documented tests/testdai
 * Renamed "Clamped BP" into "Conditioned BP"
 * Improved documentation
 * Documented tests/testdai
index c26f14a..447721a 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -101,7 +101,7 @@ ifdef WITH_CBP
 endif
 
 # Define standard libDAI header dependencies
 endif
 
 # Define standard libDAI header dependencies
-HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h $(INC)/util.h
+HEADERS=$(INC)/bipgraph.h $(INC)/graph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h $(INC)/util.h
 
 # Setup final command for C++ compiler and MEX
 ifneq ($(OS),WINDOWS)
 
 # Setup final command for C++ compiler and MEX
 ifneq ($(OS),WINDOWS)
@@ -135,6 +135,9 @@ lib: $(LIB)/libdai$(LE)
 bipgraph$(OE) : $(SRC)/bipgraph.cpp $(HEADERS)
        $(CC) -c $(SRC)/bipgraph.cpp
 
 bipgraph$(OE) : $(SRC)/bipgraph.cpp $(HEADERS)
        $(CC) -c $(SRC)/bipgraph.cpp
 
+graph$(OE) : $(SRC)/graph.cpp $(HEADERS)
+       $(CC) -c $(SRC)/graph.cpp
+
 varset$(OE) : $(SRC/varset.cpp $(HEADERS)
        $(CC) -c $(SRC)/varset.cpp
 
 varset$(OE) : $(SRC/varset.cpp $(HEADERS)
        $(CC) -c $(SRC)/varset.cpp
 
@@ -289,7 +292,7 @@ utils/fginfo$(EE) : utils/fginfo.cpp $(HEADERS) $(LIB)/libdai$(LE)
 # LIBRARY
 ##########
 
 # LIBRARY
 ##########
 
-OBJECTS:=bipgraph$(OE) varset$(OE) daialg$(OE) alldai$(OE) clustergraph$(OE) factor$(OE) factorgraph$(OE) properties$(OE) regiongraph$(OE) util$(OE) weightedgraph$(OE) exceptions$(OE) $(OBJECTS) 
+OBJECTS:=bipgraph$(OE) graph$(OE) varset$(OE) daialg$(OE) alldai$(OE) clustergraph$(OE) factor$(OE) factorgraph$(OE) properties$(OE) regiongraph$(OE) util$(OE) weightedgraph$(OE) exceptions$(OE) $(OBJECTS) 
 
 ifneq ($(OS),WINDOWS)
 $(LIB)/libdai$(LE) : $(OBJECTS)
 
 ifneq ($(OS),WINDOWS)
 $(LIB)/libdai$(LE) : $(OBJECTS)
index 7d6383a..6d438f9 100644 (file)
@@ -11,8 +11,6 @@
 /** \file
  *  \brief Contains additional doxygen documentation
  *
 /** \file
  *  \brief Contains additional doxygen documentation
  *
- *  \todo Document utils/createfg (and clean up the code a bit)
- *
  *  \idea Adapt (part of the) guidelines in http://www.boost.org/development/requirements.html#Design_and_Programming
  *
  *  \idea Use "gcc -MM" to generate dependencies for targets: http://make.paulandlesley.org/autodep.html
  *  \idea Adapt (part of the) guidelines in http://www.boost.org/development/requirements.html#Design_and_Programming
  *
  *  \idea Use "gcc -MM" to generate dependencies for targets: http://make.paulandlesley.org/autodep.html
index 50cc961..7c30825 100644 (file)
@@ -595,7 +595,7 @@ typedef TFactor<Real> Factor;
 /** \param x Variable (should be binary)
  *  \param h Field strength
  */
 /** \param x Variable (should be binary)
  *  \param h Field strength
  */
-Factor BinaryFactor( const Var &x, Real h );
+Factor createFactorIsing( const Var &x, Real h );
 
 
 /// Returns a binary pairwise factor \f$ \exp(J x_1 x_2) \f$ where \f$ x_1, x_2 = \pm 1 \f$
 
 
 /// Returns a binary pairwise factor \f$ \exp(J x_1 x_2) \f$ where \f$ x_1, x_2 = \pm 1 \f$
@@ -603,7 +603,7 @@ Factor BinaryFactor( const Var &x, Real h );
  *  \param x2 Second variable (should be binary)
  *  \param J Coupling strength
  */
  *  \param x2 Second variable (should be binary)
  *  \param J Coupling strength
  */
-Factor BinaryFactor( const Var &x1, const Var &x2, Real J );
+Factor createFactorIsing( const Var &x1, const Var &x2, Real J );
 
 
 /// Returns a random factor on the variables \a vs with strength \a beta
 
 
 /// Returns a random factor on the variables \a vs with strength \a beta
@@ -612,7 +612,7 @@ Factor BinaryFactor( const Var &x1, const Var &x2, Real J );
  *  \param vs Variables
  *  \param beta Factor strength (inverse temperature)
  */
  *  \param vs Variables
  *  \param beta Factor strength (inverse temperature)
  */
-Factor RandomFactor( const VarSet &vs, Real beta );
+Factor createFactorExpGauss( const VarSet &vs, Real beta );
 
 
 /// Returns a pairwise Potts factor \f$ \exp( J \delta_{x_1, x_2} ) \f$
 
 
 /// Returns a pairwise Potts factor \f$ \exp( J \delta_{x_1, x_2} ) \f$
@@ -620,14 +620,14 @@ Factor RandomFactor( const VarSet &vs, Real beta );
  *  \param x2 Second variable (should have the same number of states as \a x1)
  *  \param J  Factor strength
  */
  *  \param x2 Second variable (should have the same number of states as \a x1)
  *  \param J  Factor strength
  */
-Factor PottsFactor( const Var &x1, const Var &x2, Real J );
+Factor createFactorPotts( const Var &x1, const Var &x2, Real J );
 
 
 /// Returns a Kronecker delta point mass
 /** \param v Variable
  *  \param state The state of \a v that should get value 1
  */
 
 
 /// Returns a Kronecker delta point mass
 /** \param v Variable
  *  \param state The state of \a v that should get value 1
  */
-Factor DeltaFactor( const Var &v, size_t state );
+Factor createFactorDelta( const Var &v, size_t state );
 
 
 } // end of namespace dai
 
 
 } // end of namespace dai
diff --git a/include/dai/graph.h b/include/dai/graph.h
new file mode 100644 (file)
index 0000000..fe526f6
--- /dev/null
@@ -0,0 +1,288 @@
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  libDAI is licensed under the terms of the GNU General Public License version
+ *  2, or (at your option) any later version. libDAI is distributed without any
+ *  warranty. See the file COPYING for more details.
+ *
+ *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ */
+
+
+/// \file
+/// \brief Defines the GraphAL class, which represents an undirected graph as an adjacency list
+
+
+#ifndef __defined_libdai_graph_h
+#define __defined_libdai_graph_h
+
+
+#include <ostream>
+#include <vector>
+#include <algorithm>
+#include <dai/util.h>
+#include <dai/exceptions.h>
+
+
+namespace dai {
+
+
+/// Represents the neighborhood structure of nodes in an undirected graph.
+/** A graph has nodes connected by edges. Nodes are indexed by an unsigned integer. 
+ *  If there are nrNodes() nodes, they are numbered 0,1,2,...,nrNodes()-1. An edge
+ *  between node \a n1 and node \a n2 is represented by a GraphAL::Edge(\a n1,\a n2).
+ *
+ *  GraphAL is implemented as a sparse adjacency list, i.e., it stores for each node a list of
+ *  its neighboring nodes. The list of neighboring nodes is implemented as a vector of Neighbor
+ *  structures (accessible by the nb() method). Thus, each node has an associated variable of 
+ *  type GraphAL::Neighbors, which is a vector of Neighbor structures, describing its 
+ *  neighboring nodes.
+ */
+class GraphAL {
+    public:
+        /// Describes the neighbor relationship of two nodes in a GraphAL.
+        /** Sometimes we want to do an action, such as sending a
+         *  message, for all edges in a graph. However, most graphs
+         *  will be sparse, so we need some way of storing a set of
+         *  the neighbors of a node, which is both fast and
+         *  memory-efficient. We also need to be able to go between
+         *  viewing node \a a as a neighbor of node \a b, and node \a b
+         *  as a neighbor of node \a a. The Neighbor struct solves
+         *  both of these problems. Each node has a list of neighbors,
+         *  stored as a std::vector<\link Neighbor \endlink>, and 
+         *  extra information is included in the Neighbor struct which 
+         *  allows us to access a node as a neighbor of its neighbor 
+         *  (the \c dual member).
+         *
+         *  By convention, variable identifiers naming indices into a
+         *  vector of neighbors are prefixed with an underscore ("_").
+         *  The neighbor list which they point into is then understood
+         *  from context. For example:
+         *
+         *  \code
+         *  Real MR::T( size_t i, size_t _j );
+         *  \endcode
+         *
+         *  Here, \a i is the "absolute" index of node i, but \a _j is
+         *  understood as a "relative" index, giving node \a j 's entry in
+         *  <tt>nb(i)</tt>. The corresponding Neighbor structure can be
+         *  accessed as <tt>nb(i,_j)</tt> or <tt>nb(i)[_j]</tt>. The 
+         *  absolute index of \a _j, which would be called \a j, can be 
+         *  recovered from the \c node member: <tt>nb(i,_j).node</tt>. 
+         *  The \c iter member gives the relative index \a _j, and the 
+         *  \c dual member gives the "dual" relative index, i.e., the 
+         *  index of \a i in \a j 's neighbor list.
+         *
+         *  \code
+         *  Neighbor n = nb(i,_j);
+         *  n.node == j &&
+         *  n.iter == _j &&
+         *  nb(n.node,n.dual).node == i
+         *  \endcode
+         *
+         *  There is no easy way to transform a pair of absolute node
+         *  indices \a i and \a j into a Neighbor structure relative
+         *  to one of the nodes. Such a feature has never yet been
+         *  found to be necessary. Iteration over edges can always be
+         *  accomplished using the Neighbor lists, and by writing
+         *  functions that accept relative indices:
+         *  \code
+         *  for( size_t i = 0; i < nrNodes(); ++i )
+         *      foreach( const Neighbor &j, nb(i) )
+         *          T( i, j.iter );
+         *  \endcode
+         */
+        struct Neighbor {
+            /// Corresponds to the index of this Neighbor entry in the vector of neighbors
+            size_t iter;
+            /// Contains the number of the neighboring node
+            size_t node;
+            /// Contains the "dual" iter
+            size_t dual;
+
+            /// Default constructor
+            Neighbor() {}
+            /// Constructor that sets the Neighbor members according to the parameters
+            Neighbor( size_t iter, size_t node, size_t dual ) : iter(iter), node(node), dual(dual) {}
+
+            /// Cast to \c size_t returns \c node member
+            operator size_t () const { return node; }
+        };
+
+        /// Describes the neighbors of some node.
+        typedef std::vector<Neighbor> Neighbors;
+
+        /// Represents an edge: an Edge(\a n1,\a n2) corresponds to the edge between node \a n1 and node \a n2.
+        typedef std::pair<size_t,size_t> Edge;
+
+    private:
+        /// Contains for each node a vector of its neighbors
+        std::vector<Neighbors> _nb;
+
+        /// Used internally by isTree()
+        typedef std::vector<size_t> levelType;
+
+    public:
+    /// \name Constructors and destructors
+    //@{
+        /// Default constructor (creates an empty graph).
+        GraphAL() : _nb() {}
+
+        /// Constructs GraphAL with \a nr nodes and no edges.
+        GraphAL( size_t nr ) : _nb( nr ) {}
+
+        /// Constructs GraphAL from a range of edges.
+        /** \tparam EdgeInputIterator Iterator that iterates over instances of GraphAL::Edge.
+         *  \param nr The number of nodes.
+         *  \param begin Points to the first edge.
+         *  \param end Points just beyond the last edge.
+         */
+        template<typename EdgeInputIterator>
+        GraphAL( size_t nr, EdgeInputIterator begin, EdgeInputIterator end ) : _nb() {
+            construct( nr, begin, end );
+        }
+    //@}
+
+    /// \name Accessors and mutators
+    //@{
+        /// Returns constant reference to the \a _n2 'th neighbor of node \a n1
+        const Neighbor & nb( size_t n1, size_t _n2 ) const {
+            DAI_DEBASSERT( n1 < _nb.size() );
+            DAI_DEBASSERT( _n2 < _nb[n1].size() );
+            return _nb[n1][_n2];
+        }
+        /// Returns reference to the \a _n2 'th neighbor of node \a n1
+        Neighbor & nb( size_t n1, size_t _n2 ) {
+            DAI_DEBASSERT( n1 < _nb.size() );
+            DAI_DEBASSERT( _n2 < _nb[n1].size() );
+            return _nb[n1][_n2];
+        }
+
+        /// Returns constant reference to all neighbors of node \a n
+        const Neighbors & nb( size_t n ) const {
+            DAI_DEBASSERT( n < _nb.size() );
+            return _nb[n];
+        }
+        /// Returns reference to all neighbors of node \a n
+        Neighbors & nb( size_t n ) {
+            DAI_DEBASSERT( n < _nb.size() );
+            return _nb[n];
+        }
+    //@}
+
+    /// \name Adding nodes and edges
+    //@{
+        /// (Re)constructs GraphAL from a range of edges.
+        /** \tparam EdgeInputIterator Iterator that iterates over instances of GraphAL::Edge.
+         *  \param nr The number of nodes.
+         *  \param begin Points to the first edge.
+         *  \param end Points just beyond the last edge.
+         */
+        template<typename EdgeInputIterator>
+        void construct( size_t nr, EdgeInputIterator begin, EdgeInputIterator end );
+
+        /// Adds a node without neighbors and returns the index of the added node.
+        size_t addNode() { _nb.push_back( Neighbors() ); return _nb.size() - 1; }
+
+        /// Adds a node, with neighbors specified by a range of nodes.
+        /** \tparam NodeInputIterator Iterator that iterates over instances of \c size_t.
+         *  \param begin Points to the first index of the nodes that should become neighbors of the added node.
+         *  \param end Points just beyond the last index of the nodes that should become neighbors of the added node.
+         *  \param sizeHint For improved efficiency, the size of the range may be specified by \a sizeHint.
+         *  \returns Index of the added node.
+         */
+        template <typename NodeInputIterator>
+        size_t addNode( NodeInputIterator begin, NodeInputIterator end, size_t sizeHint = 0 ) {
+            Neighbors nbsnew;
+            nbsnew.reserve( sizeHint );
+            size_t iter = 0;
+            for( NodeInputIterator it = begin; it != end; ++it ) {
+                DAI_ASSERT( *it < nrNodes() );
+                Neighbor nb1new( iter, *it, nb(*it).size() );
+                Neighbor nb2new( nb(*it).size(), nrNodes(), iter++ );
+                nbsnew.push_back( nb1new );
+                nb( *it ).push_back( nb2new );
+            }
+            _nb.push_back( nbsnew );
+            return _nb.size() - 1;
+        }
+
+        /// Adds an edge between node \a n1 and node \a n2.
+        /** If \a check == \c true, only adds the edge if it does not exist already.
+         */
+        void addEdge( size_t n1, size_t n2, bool check = true );
+    //@}
+
+    /// \name Erasing nodes and edges
+    //@{
+        /// Removes node \a n and all incident edges; indices of other nodes are changed accordingly.
+        void eraseNode( size_t n );
+
+        /// Removes edge between node \a n1 and node \a n2.
+        void eraseEdge( size_t n1, size_t n2 );
+    //@}
+
+    /// \name Queries
+    //@{
+        /// Returns number of nodes
+        size_t nrNodes() const { return _nb.size(); }
+
+        /// Calculates the number of edges, time complexity: O(nrNodes())
+        size_t nrEdges() const {
+            size_t sum = 0;
+            for( size_t i = 0; i < nrNodes(); i++ )
+                sum += nb(i).size();
+            return sum;
+        }
+
+        /// Returns true if the graph is connected
+        bool isConnected() const;
+
+        /// Returns true if the graph is a tree, i.e., if it is singly connected and connected.
+        bool isTree() const;
+
+        /// Checks internal consistency
+        void checkConsistency() const;
+    //@}
+
+    /// \name Input and output
+    //@{
+        /// Writes this GraphAL to an output stream in GraphALViz .dot syntax
+        void printDot( std::ostream& os ) const;
+    //@}
+};
+
+
+template<typename EdgeInputIterator>
+void GraphAL::construct( size_t nr, EdgeInputIterator begin, EdgeInputIterator end ) {
+    _nb.clear();
+    _nb.resize( nr );
+
+    for( EdgeInputIterator e = begin; e != end; e++ ) {
+#ifdef DAI_DEBUG
+        addEdge( e->first, e->second, true );
+#else
+        addEdge( e->first, e->second, false );
+#endif
+    }
+}
+
+
+/// Creates a fully-connected graph with \a N nodes
+GraphAL createGraphFull( size_t N );
+/// Creates a two-dimensional rectangular grid of \a n1 by \a n2 nodes, which can be \a periodic
+GraphAL createGraphGrid( size_t n1, size_t n2, bool periodic );
+/// Creates a three-dimensional rectangular grid of \a n1 by \a n2 by \a n3 nodes, which can be \a periodic
+GraphAL createGraphGrid3D( size_t n1, size_t n2, size_t n3, bool periodic );
+/// Creates a graph consisting of a single loop of \a N nodes
+GraphAL createGraphLoop( size_t N );
+/// Creates a random tree-structured graph of \a N nodes
+GraphAL createGraphTree( size_t N );
+/// Creates a random regular graph of \a N nodes with uniform connectivity \a d
+GraphAL createGraphRegular( size_t N, size_t d );
+
+
+} // end of namespace dai
+
+
+#endif
index cc723bc..0d543b5 100644 (file)
@@ -71,9 +71,15 @@ class DEdge {
 class UEdge {
     public:
         /// First node index
 class UEdge {
     public:
         /// First node index
-        size_t  n1;
+        union {
+            size_t n1;
+            size_t first;   /// alias
+        };
         /// Second node index
         /// Second node index
-        size_t  n2;
+        union {
+            size_t n2;
+            size_t second;  /// alias
+        };
 
         /// Default constructor
         UEdge() {}
 
         /// Default constructor
         UEdge() {}
@@ -112,19 +118,23 @@ class UEdge {
 
 
 /// Represents an undirected graph, implemented as a std::set of undirected edges
 
 
 /// Represents an undirected graph, implemented as a std::set of undirected edges
-class Graph : public std::set<UEdge> {
+class GraphEL : public std::set<UEdge> {
     public:
         /// Default constructor
     public:
         /// Default constructor
-        Graph() {}
+        GraphEL() {}
 
         /// Construct from range of objects that can be cast to DEdge
         template <class InputIterator>
 
         /// Construct from range of objects that can be cast to DEdge
         template <class InputIterator>
-        Graph( InputIterator begin, InputIterator end ) {
+        GraphEL( InputIterator begin, InputIterator end ) {
             insert( begin, end );
         }
 };
 
 
             insert( begin, end );
         }
 };
 
 
+/// \deprecated Please use GraphEL instead.
+typedef GraphEL Graph;
+
+
 /// Represents an undirected weighted graph, with weights of type \a T, implemented as a std::map mapping undirected edges to weights
 template<class T> class WeightedGraph : public std::map<UEdge, T> {};
 
 /// Represents an undirected weighted graph, with weights of type \a T, implemented as a std::map mapping undirected edges to weights
 template<class T> class WeightedGraph : public std::map<UEdge, T> {};
 
@@ -142,7 +152,7 @@ class RootedTree : public std::vector<DEdge> {
         /// Constructs a rooted tree from a tree and a root
         /** \pre T has no cycles and contains node \a Root
          */
         /// Constructs a rooted tree from a tree and a root
         /** \pre T has no cycles and contains node \a Root
          */
-        RootedTree( const Graph &T, size_t Root );
+        RootedTree( const GraphEL &T, size_t Root );
 };
 
 
 };
 
 
@@ -174,7 +184,7 @@ template<typename T> RootedTree MinSpanningTreePrims( const WeightedGraph<T> &G
         prim_minimum_spanning_tree( g, &(p[0]) );
 
         // Store tree edges in result
         prim_minimum_spanning_tree( g, &(p[0]) );
 
         // Store tree edges in result
-        Graph tree;
+        GraphEL tree;
         size_t root = 0;
         for( size_t i = 0; i != p.size(); i++ )
             if( p[i] != i )
         size_t root = 0;
         for( size_t i = 0; i != p.size(); i++ )
             if( p[i] != i )
@@ -217,7 +227,7 @@ template<typename T> RootedTree MaxSpanningTreePrims( const WeightedGraph<T> &G
  *  (which becomes uniform in the limit that \a d is small and \a N goes
  *  to infinity).
  */
  *  (which becomes uniform in the limit that \a d is small and \a N goes
  *  to infinity).
  */
-Graph RandomDRegularGraph( size_t N, size_t d );
+GraphEL RandomDRegularGraph( size_t N, size_t d );
 
 
 } // end of namespace dai
 
 
 } // end of namespace dai
index b4f69de..46c73fa 100644 (file)
@@ -18,7 +18,7 @@ namespace dai {
 using namespace std;
 
 
 using namespace std;
 
 
-Factor BinaryFactor( const Var &n, Real h ) {
+Factor createFactorIsing( const Var &n, Real h ) {
     DAI_ASSERT( n.states() == 2 );
     Real buf[2];
     buf[0] = std::exp(-h);
     DAI_ASSERT( n.states() == 2 );
     Real buf[2];
     buf[0] = std::exp(-h);
@@ -27,7 +27,7 @@ Factor BinaryFactor( const Var &n, Real h ) {
 }
 
 
 }
 
 
-Factor BinaryFactor( const Var &n1, const Var &n2, Real J ) {
+Factor createFactorIsing( const Var &n1, const Var &n2, Real J ) {
     DAI_ASSERT( n1.states() == 2 );
     DAI_ASSERT( n2.states() == 2 );
     DAI_ASSERT( n1 != n2 );
     DAI_ASSERT( n1.states() == 2 );
     DAI_ASSERT( n2.states() == 2 );
     DAI_ASSERT( n1 != n2 );
@@ -38,7 +38,7 @@ Factor BinaryFactor( const Var &n1, const Var &n2, Real J ) {
 }
 
 
 }
 
 
-Factor RandomFactor( const VarSet &ns, Real beta ) {
+Factor createFactorExpGauss( const VarSet &ns, Real beta ) {
     Factor fac( ns );
     for( size_t t = 0; t < fac.states(); t++ )
         fac[t] = std::exp(rnd_stdnormal() * beta);
     Factor fac( ns );
     for( size_t t = 0; t < fac.states(); t++ )
         fac[t] = std::exp(rnd_stdnormal() * beta);
@@ -46,7 +46,7 @@ Factor RandomFactor( const VarSet &ns, Real beta ) {
 }
 
 
 }
 
 
-Factor PottsFactor( const Var &n1, const Var &n2, Real J ) {
+Factor createFactorPotts( const Var &n1, const Var &n2, Real J ) {
     Factor fac( VarSet( n1, n2 ), 1.0 );
     DAI_ASSERT( n1.states() == n2.states() );
     for( size_t s = 0; s < n1.states(); s++ )
     Factor fac( VarSet( n1, n2 ), 1.0 );
     DAI_ASSERT( n1.states() == n2.states() );
     for( size_t s = 0; s < n1.states(); s++ )
@@ -55,7 +55,7 @@ Factor PottsFactor( const Var &n1, const Var &n2, Real J ) {
 }
 
 
 }
 
 
-Factor DeltaFactor( const Var &v, size_t state ) {
+Factor createFactorDelta( const Var &v, size_t state ) {
     Factor fac( v, 0.0 );
     DAI_ASSERT( state < v.states() );
     fac[state] = 1.0;
     Factor fac( v, 0.0 );
     DAI_ASSERT( state < v.states() );
     fac[state] = 1.0;
diff --git a/src/graph.cpp b/src/graph.cpp
new file mode 100644 (file)
index 0000000..73e5a54
--- /dev/null
@@ -0,0 +1,306 @@
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  libDAI is licensed under the terms of the GNU General Public License version
+ *  2, or (at your option) any later version. libDAI is distributed without any
+ *  warranty. See the file COPYING for more details.
+ *
+ *  Copyright (C) 2006-2010  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ */
+
+
+#include <dai/graph.h>
+#include <dai/weightedgraph.h>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/connected_components.hpp>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+void GraphAL::addEdge( size_t n1, size_t n2, bool check ) {
+    DAI_ASSERT( n1 < nrNodes() );
+    DAI_ASSERT( n2 < nrNodes() );
+    bool exists = false;
+    if( check ) {
+        // Check whether the edge already exists
+        foreach( const Neighbor &n, nb(n1) )
+            if( n == n2 ) {
+                exists = true;
+                break;
+            }
+    }
+    if( !exists ) { // Add edge
+        Neighbor nb_1( nb(n1).size(), n2, nb(n2).size() );
+        Neighbor nb_2( nb_1.dual, n1, nb_1.iter );
+        nb(n1).push_back( nb_1 );
+        nb(n2).push_back( nb_2 );
+    }
+}
+
+
+void GraphAL::eraseNode( size_t n ) {
+    DAI_ASSERT( n < nrNodes() );
+    // Erase neighbor entry of node n
+    _nb.erase( _nb.begin() + n );
+    // Adjust neighbor entries of nodes
+    for( size_t n2 = 0; n2 < nrNodes(); n2++ ) {
+        for( size_t iter = 0; iter < nb(n2).size(); ) {
+            Neighbor &m = nb(n2, iter);
+            if( m.node == n ) {
+                // delete this entry, because it points to the deleted node
+                nb(n2).erase( nb(n2).begin() + iter );
+            } else if( m.node > n ) {
+                // update this entry and the corresponding dual of the neighboring node
+                m.iter = iter;
+                m.node--;
+                nb( m.node, m.dual ).dual = iter;
+                iter++;
+            } else {
+                // skip
+                iter++;
+            }
+        }
+    }
+}
+
+
+void GraphAL::eraseEdge( size_t n1, size_t n2 ) {
+    DAI_ASSERT( n1 < nrNodes() );
+    DAI_ASSERT( n2 < nrNodes() );
+    size_t iter;
+    // Search for edge among neighbors of n1
+    for( iter = 0; iter < nb(n1).size(); iter++ )
+        if( nb(n1, iter).node == n2 ) {
+            // Remove it
+            nb(n1).erase( nb(n1).begin() + iter );
+            break;
+        }
+    // Change the iter and dual values of the subsequent neighbors
+    for( ; iter < nb(n1).size(); iter++ ) {
+        Neighbor &m = nb( n1, iter );
+        m.iter = iter;
+        nb( m.node, m.dual ).dual = iter;
+    }
+    // Search for edge among neighbors of n2
+    for( iter = 0; iter < nb(n2).size(); iter++ )
+        if( nb(n2, iter).node == n1 ) {
+            // Remove it
+            nb(n2).erase( nb(n2).begin() + iter );
+            break;
+        }
+    // Change the iter and node values of the subsequent neighbors
+    for( ; iter < nb(n2).size(); iter++ ) {
+        Neighbor &m = nb( n2, iter );
+        m.iter = iter;
+        nb( m.node, m.dual ).dual = iter;
+    }
+}
+
+
+bool GraphAL::isConnected() const {
+    if( nrNodes() == 0 ) {
+        return true;
+    } else {
+        std::vector<bool> incomponent( nrNodes(), false );
+
+        incomponent[0] = true;
+        bool found_new_nodes;
+        do {
+            found_new_nodes = false;
+
+            // For all nodes, check if they are connected with the (growing) component
+            for( size_t n1 = 0; n1 < nrNodes(); n1++ )
+                if( !incomponent[n1] ) {
+                    foreach( const Neighbor &n2, nb(n1) ) {
+                        if( incomponent[n2] ) {
+                            found_new_nodes = true;
+                            incomponent[n1] = true;
+                            break;
+                        }
+                    }
+                }
+        } while( found_new_nodes );
+
+        // Check if there are remaining nodes (not in the component)
+        bool all_connected = true;
+        for( size_t n1 = 0; (n1 < nrNodes()) && all_connected; n1++ )
+            if( !incomponent[n1] )
+                all_connected = false;
+
+        return all_connected;
+
+        // BGL implementation is slower...
+    /*  using namespace boost;
+        typedef adjacency_list< vecS, vecS, undirectedS, property<vertex_distance_t, int> > boostGraphAL;
+        typedef pair<size_t, size_t> E;
+
+        // Copy graph structure into boostGraphAL object
+        vector<E> edges;
+        edges.reserve( nrEdges() );
+        for( size_t n1 = 0; n1 < nrNodes(); n1++ )
+            foreach( const Neighbor &n2, nb(n1) )
+                if( n1 < n2 )
+                    edges.push_back( E( n1, n2 ) );
+        boostGraphAL g( edges.begin(), edges.end(), nrNodes() );
+
+        // Construct connected components using Boost GraphAL Library
+        std::vector<int> component( num_vertices( g ) );
+        int num_comp = connected_components( g, make_iterator_property_map(component.begin(), get(vertex_index, g)) );
+
+        return (num_comp == 1);
+    */
+    }
+}
+
+
+bool GraphAL::isTree() const {
+    using namespace std;
+    vector<levelType> levels;
+
+    bool foundCycle = false;
+    size_t Nr = 0;
+
+    if( nrNodes() == 0 )
+        return true;
+    else {
+        levelType newLevel;
+        do {
+            newLevel.clear();
+            if( levels.size() == 0 ) {
+                size_t n1 = 0;
+                // add n1
+                newLevel = vector<size_t>( 1, n1 );
+            } else {
+                const levelType &prevLevel = levels.back();
+                // build newLevel
+                foreach( size_t n2, prevLevel ) { // for all n2 in the previous level
+                    foreach( const Neighbor &n1, nb(n2) ) { // for all neighbors n1 of n2
+                        if( find( prevLevel.begin(), prevLevel.end(), n1 ) == prevLevel.end() ) { // n1 not in previous level
+                            if( find( newLevel.begin(), newLevel.end(), n1 ) != newLevel.end() )
+                                foundCycle = true; // n1 already in new level: we found a cycle
+                            else
+                                newLevel.push_back( n1 ); // add n1 to new level
+                        }
+                        if( foundCycle )
+                            break;
+                    }
+                    if( foundCycle )
+                        break;
+                }
+            }
+            levels.push_back( newLevel );
+            Nr += newLevel.size();
+        } while( (newLevel.size() != 0) && !foundCycle );
+        if( Nr == nrNodes() && !foundCycle )
+            return true;
+        else
+            return false;
+    }
+}
+
+
+void GraphAL::printDot( std::ostream& os ) const {
+    using namespace std;
+    os << "graph G {" << endl;
+    os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
+    for( size_t n = 0; n < nrNodes(); n++ )
+        os << "\tx" << n << ";" << endl;
+    for( size_t n1 = 0; n1 < nrNodes(); n1++ )
+        foreach( const Neighbor &n2, nb(n1) )
+            os << "\tx" << n1 << " -- y" << n2 << ";" << endl;
+    os << "}" << endl;
+}
+
+
+void GraphAL::checkConsistency() const {
+    size_t N = nrNodes();
+    for( size_t n1 = 0; n1 < N; n1++ ) {
+        size_t iter = 0;
+        foreach( const Neighbor &n2, nb(n1) ) {
+            DAI_ASSERT( n2.iter == iter );
+            DAI_ASSERT( n2.node < N );
+            DAI_ASSERT( n2.dual < nb(n2).size() );
+            DAI_ASSERT( nb(n2, n2.dual) == n1 );
+            iter++;
+        }
+    }
+    for( size_t n2 = 0; n2 < N; n2++ ) {
+        size_t iter = 0;
+        foreach( const Neighbor &n1, nb(n2) ) {
+            DAI_ASSERT( n1.iter == iter );
+            DAI_ASSERT( n1.node < N );
+            DAI_ASSERT( n1.dual < nb(n1).size() );
+            DAI_ASSERT( nb(n1, n1.dual) == n2 );
+            iter++;
+        }
+    }
+}
+
+
+GraphAL createGraphFull( size_t N ) {
+    GraphAL result( N );
+    for( size_t i = 0; i < N; i++ )
+        for( size_t j = i+1; j < N; j++ )
+            result.addEdge( i, j, false );
+    return result;
+}
+
+
+GraphAL createGraphGrid( size_t n1, size_t n2, bool periodic ) {
+    GraphAL result( n1*n2 );
+    for( size_t i1 = 0; i1 < n1; i1++ )
+        for( size_t i2 = 0; i2 < n2; i2++ ) {
+            if( i1+1 < n1 || periodic )
+                result.addEdge( i1*n2 + i2, ((i1+1)%n1)*n2 + i2, n1 <= 2 );
+            if( i2+1 < n2 || periodic )
+                result.addEdge( i1*n2 + i2, i1*n2 + ((i2+1)%n2), n2 <= 2 );
+        }
+    return result;
+}
+
+
+GraphAL createGraphGrid3D( size_t n1, size_t n2, size_t n3, bool periodic ) {
+    GraphAL result( n1*n2*n3 );
+    for( size_t i1 = 0; i1 < n1; i1++ )
+        for( size_t i2 = 0; i2 < n2; i2++ )
+            for( size_t i3 = 0; i3 < n3; i3++ ) {
+                if( i1+1 < n1 || periodic )
+                    result.addEdge( i1*n2*n3 + i2*n3 + i3, ((i1+1)%n1)*n2*n3 + i2*n3 + i3, n1 <= 2 );
+                if( i2+1 < n2 || periodic )
+                    result.addEdge( i1*n2*n3 + i2*n3 + i3, i1*n2*n3 + ((i2+1)%n2)*n3 + i3, n2 <= 2 );
+                if( i3+1 < n2 || periodic )
+                    result.addEdge( i1*n2*n3 + i2*n3 + i3, i1*n2*n3 + i2*n3 + ((i3+1)%n3), n3 <= 2 );
+            }
+    return result;
+}
+
+
+GraphAL createGraphLoop( size_t N ) {
+    GraphAL result( N );
+    for( size_t i = 0; i < N; i++ )
+        result.addEdge( i, (i+1)%N, N <= 2 );
+    return result;
+}
+
+
+GraphAL createGraphTree( size_t N ) {
+    GraphAL result( N );
+    for( size_t i = 1; i < N; i++ ) {
+        size_t j = rnd_int( 0, i-1 );
+        result.addEdge( i, j, false );
+    }
+    return result;
+}
+
+
+GraphAL createGraphRegular( size_t N, size_t d ) {
+    GraphEL g = RandomDRegularGraph( N, d );
+    return GraphAL( N, g.begin(), g.end() );
+}
+
+
+} // end of namespace dai
index 31a448e..b1d5edd 100644 (file)
@@ -381,7 +381,7 @@ size_t JTree::findEfficientTree( const VarSet& vs, RootedTree &Tree, size_t Prev
     }
 
     // reorder the tree edges such that maxalpha becomes the new root
     }
 
     // reorder the tree edges such that maxalpha becomes the new root
-    RootedTree newTree( Graph( RTree.begin(), RTree.end() ), maxalpha );
+    RootedTree newTree( GraphEL( RTree.begin(), RTree.end() ), maxalpha );
 
     // identify subtree that contains all variables of vs which are not in the new root
     set<DEdge> subTree;
 
     // identify subtree that contains all variables of vs which are not in the new root
     set<DEdge> subTree;
index 6f199c7..e70e4bf 100644 (file)
@@ -21,10 +21,10 @@ namespace dai {
 using namespace std;
 
 
 using namespace std;
 
 
-RootedTree::RootedTree( const Graph &T, size_t Root ) {
+RootedTree::RootedTree( const GraphEL &T, size_t Root ) {
     if( T.size() != 0 ) {
         // Make a copy
     if( T.size() != 0 ) {
         // Make a copy
-        Graph Gr = T;
+        GraphEL Gr = T;
 
         // Nodes in the tree
         set<size_t> nodes;
 
         // Nodes in the tree
         set<size_t> nodes;
@@ -34,7 +34,7 @@ RootedTree::RootedTree( const Graph &T, size_t Root ) {
 
         // Keep adding edges until done
         while( !(Gr.empty()) )
 
         // Keep adding edges until done
         while( !(Gr.empty()) )
-            for( Graph::iterator e = Gr.begin(); e != Gr.end(); ) {
+            for( GraphEL::iterator e = Gr.begin(); e != Gr.end(); ) {
                 bool e1_in_nodes = nodes.count( e->n1 );
                 bool e2_in_nodes = nodes.count( e->n2 );
                 DAI_ASSERT( !(e1_in_nodes && e2_in_nodes) );
                 bool e1_in_nodes = nodes.count( e->n1 );
                 bool e2_in_nodes = nodes.count( e->n2 );
                 DAI_ASSERT( !(e1_in_nodes && e2_in_nodes) );
@@ -57,7 +57,7 @@ RootedTree::RootedTree( const Graph &T, size_t Root ) {
 }
 
 
 }
 
 
-Graph RandomDRegularGraph( size_t N, size_t d ) {
+GraphEL RandomDRegularGraph( size_t N, size_t d ) {
     DAI_ASSERT( (N * d) % 2 == 0 );
 
     bool ready = false;
     DAI_ASSERT( (N * d) % 2 == 0 );
 
     bool ready = false;
@@ -117,7 +117,7 @@ Graph RandomDRegularGraph( size_t N, size_t d ) {
             ready = false;
     }
 
             ready = false;
     }
 
-    return Graph( G.begin(), G.end() );
+    return GraphEL( G.begin(), G.end() );
 }
 
 
 }
 
 
index 9dfd1ca..1685526 100644 (file)
@@ -4,7 +4,7 @@
  *  2, or (at your option) any later version. libDAI is distributed without any
  *  warranty. See the file COPYING for more details.
  *
  *  2, or (at your option) any later version. libDAI is distributed without any
  *  warranty. See the file COPYING for more details.
  *
- *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2006-2010  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
  */
 
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
  */
 
 #include <iterator>
 #include <algorithm>
 #include <boost/program_options.hpp>
 #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>
 #include <dai/factorgraph.h>
 #include <dai/weightedgraph.h>
 #include <dai/util.h>
+#include <dai/graph.h>
+#include <dai/enum.h>
+#include <dai/properties.h>
 
 
 using namespace std;
 using namespace dai;
 namespace po = boost::program_options;
 
 
 using namespace std;
 using namespace dai;
 namespace po = boost::program_options;
-typedef boost::numeric::ublas::compressed_matrix<Real> matrix;
-typedef matrix::value_array_type::const_iterator matrix_vcit;
-typedef matrix::index_array_type::const_iterator matrix_icit;
 
 
 
 
-// w should be upper triangular or lower triangular
-void WTh2FG( const matrix &w, const vector<Real> &th, FactorGraph &fg ) {
-    vector<Var>    vars;
-    vector<Factor> factors;
+/// Possible factor types
+DAI_ENUM(FactorType,ISING,EXPGAUSS,POTTS);
 
 
-    size_t N = th.size();
-    DAI_ASSERT( (w.size1() == N) && (w.size2() == N) );
 
 
-    vars.reserve(N);
+/// Creates a factor graph from a pairwise interactions graph
+/** \param G  Graph describing interactions between variables
+ *  \param ft Type of factors to use for interactions
+ *  \param states Number of states of the variables
+ *  \param props Additional properties for generating the interactions
+ */
+FactorGraph createFG( const GraphAL &G, FactorType ft, size_t states, const PropertySet &props ) {
+    size_t N = G.nrNodes();
+
+    if( states > 2 )
+        DAI_ASSERT( ft != FactorType::ISING );
+
+    // Get inverse temperature
+    Real beta = 1.0;
+    if( ft != FactorType::ISING ) 
+        beta = props.GetAs<Real>("beta");
+
+    // Get properties for Ising factors
+    Real mean_h = 0.0;
+    Real sigma_h = 0.0;
+    Real mean_J = 0.0;
+    Real sigma_J = 0.0;
+    if( ft == FactorType::ISING ) {
+        mean_h = props.GetAs<Real>("mean_th");
+        sigma_h = props.GetAs<Real>("sigma_th");
+        mean_J = props.GetAs<Real>("mean_w");
+        sigma_J = props.GetAs<Real>("sigma_w");
+    }
+    
+    // Create variables
+    vector<Var> vars;
+    vars.reserve( N );
     for( size_t i = 0; i < N; i++ )
     for( size_t i = 0; i < N; i++ )
-        vars.push_back(Var(i,2));
+        vars.push_back( Var( i, states ) );
 
 
-    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];
-        Real w_ij = w.value_data()[pos];
-        factors.push_back( BinaryFactor( vars[i], vars[j], w_ij ) );
-    }
+    // Create factors
+    vector<Factor> factors;
+    factors.reserve( G.nrEdges() + N );
+    // Pairwise factors
     for( size_t i = 0; i < N; i++ )
     for( size_t i = 0; i < N; i++ )
-        factors.push_back( BinaryFactor( vars[i], th[i] ) );
+        foreach( const GraphAL::Neighbor &j, G.nb(i) )
+            if( i < j ) {
+                if( ft == FactorType::POTTS )
+                    factors.push_back( createFactorPotts( vars[i], vars[j], beta ) );
+                else if( ft == FactorType::EXPGAUSS )
+                    factors.push_back( createFactorExpGauss( VarSet( vars[i], vars[j] ), beta ) );
+                else if( ft == FactorType::ISING ) {
+                    Real J = rnd_stdnormal() * sigma_J + mean_J;
+                    factors.push_back( createFactorIsing( vars[i], vars[j], J ) );
+                }
+            }
+    // Unary factors
+    if( ft == FactorType::ISING )
+        for( size_t i = 0; i < N; i++ ) {
+            Real h = rnd_stdnormal() * sigma_h + mean_h;
+            factors.push_back( createFactorIsing( vars[i], h ) );
+        }
 
 
-    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+    return FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
 }
 
 
 }
 
 
-void MakeHOIFG( size_t N, size_t M, size_t k, Real sigma, FactorGraph &fg ) {
+/// Return a random factor graph with higher-order interactions
+/** \param N number of variables
+ *  \param M number of factors
+ *  \param k number of variables that each factor depends on
+ *  \param beta standard-deviation of Gaussian log-factor entries
+ */
+FactorGraph createHOIFG( size_t N, size_t M, size_t k, Real beta ) {
     vector<Var> vars;
     vector<Factor> factors;
 
     vector<Var> vars;
     vector<Factor> factors;
 
@@ -82,172 +120,36 @@ void MakeHOIFG( size_t N, size_t M, size_t k, Real sigma, FactorGraph &fg ) {
                 }
             } while( 1 );
         }
                 }
             } while( 1 );
         }
-        factors.push_back( RandomFactor( vars, sigma ) );
-    }
-
-    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
-}
-
-
-void MakeFullFG( size_t N, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) {
-    matrix w(N,N,N*(N-1)/2);
-    vector<Real> th(N,0.0);
-
-    for( size_t i = 0; i < N; i++ ) {
-        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;
-    }
-
-    WTh2FG( w, th, fg );
-}
-
-
-void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, Real 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 ) );
-            }
-
-    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
-}
-
-
-void MakeGridFG( long periodic, size_t n, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) {
-    size_t N = n*n;
-
-    matrix w(N,N,2*N);
-    vector<Real> 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) = rnd_stdnormal() * sigma_w + mean_w;
-            if( j+1 < n || periodic )
-                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;
-        }
-
-    WTh2FG( w, th, fg );
-}
-
-
-void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, Real 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( VarSet( vars[i*n+j], vars[((i+1)%n)*n+j] ), beta ) );
-            if( j+1 < n || periodic )
-                factors.push_back( RandomFactor( VarSet( 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() );
-}
-
-
-void MakeLoopFG( size_t N, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) {
-    matrix w(N,N,N);
-    vector<Real> 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;
-    }
-
-    WTh2FG( w, th, fg );
-}
-
-
-void MakeLoopNonbinaryFG( size_t N, size_t states, Real beta, FactorGraph &fg ) {
-    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( N );
-    for( size_t i = 0; i < N; i++ ) {
-        factors.push_back( RandomFactor( VarSet( vars[i], vars[(i+1)%N] ), beta ) );
+        factors.push_back( createFactorExpGauss( vars, beta ) );
     }
 
     }
 
-    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
+    return FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
 }
 
 
 }
 
 
-void MakeTreeFG( size_t N, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) {
-    matrix w(N,N,N-1);
-    vector<Real> th(N,0.0);
-
-    for( size_t i = 0; i < N; i++ ) {
-        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;
-        }
-    }
-
-    WTh2FG( w, th, fg );
-}
-
-
-void MakeDRegFG( size_t N, size_t d, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) {
-    matrix w(N,N,(d*N)/2);
-    vector<Real> th(N,0.0);
-
-    Graph g = RandomDRegularGraph( N, d );
-    foreach( const UEdge &e, g )
-        w(e.n1, e.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 ) {
+/// Creates a regular random bipartite graph
+/** \param N1 = number of nodes of type 1
+ *  \param d1 = size of neighborhoods of nodes of type 1
+ *  \param N2 = number of nodes of type 2
+ *  \param d2 = size of neighborhoods of nodes of type 2
+ *  \note asserts that N1 * d1 == N2 * d2
+ */
+BipartiteGraph createRandomBipartiteGraph( size_t N1, size_t N2, size_t d1, size_t d2 ) {
     BipartiteGraph G;
 
     BipartiteGraph G;
 
-    DAI_ASSERT( N * n == K * k );
+    DAI_ASSERT( N1 * d1 == N2 * d2 );
 
     // build lists of degree-repeated vertex numbers
 
     // 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;
+    std::vector<size_t> stubs1( N1*d1, 0 );
+    for( size_t n1 = 0; n1 < N1; n1++ )
+        for( size_t t = 0; t < d1; t++ )
+            stubs1[n1*d1 + t] = n1;
 
     // build lists of degree-repeated vertex numbers
 
     // 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;
+    std::vector<size_t> stubs2( N2*d2, 0 );
+    for( size_t n2 = 0; n2 < N2; n2++ )
+        for( size_t t = 0; t < d2; t++ )
+            stubs2[n2*d2 + t] = n2;
 
     // shuffle lists
     random_shuffle( stubs1.begin(), stubs1.end() );
 
     // shuffle lists
     random_shuffle( stubs1.begin(), stubs1.end() );
@@ -255,18 +157,18 @@ BipartiteGraph CreateRandomBipartiteGraph( size_t N, size_t K, size_t n, size_t
 
     // add edges
     vector<BipartiteGraph::Edge> edges;
 
     // add edges
     vector<BipartiteGraph::Edge> edges;
-    edges.reserve( N*n );
-    for( size_t e = 0; e < N*n; e++ )
+    edges.reserve( N1*d1 );
+    for( size_t e = 0; e < N1*d1; e++ )
         edges.push_back( BipartiteGraph::Edge(stubs1[e], stubs2[e]) );
 
     // finish construction
         edges.push_back( BipartiteGraph::Edge(stubs1[e], stubs2[e]) );
 
     // finish construction
-    G.construct( N, K, edges.begin(), edges.end() );
+    G.construct( N1, N2, edges.begin(), edges.end() );
 
     return G;
 }
 
 
 
     return G;
 }
 
 
-// Returns x**n % p, assuming p is prime
+/// 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++ )
 int powmod (int x, int n, int p) {
     int y = 1;
     for( int m = 0; m < n; m++ )
@@ -275,7 +177,7 @@ int powmod (int x, int n, int p) {
 }
 
 
 }
 
 
-// Returns order of x in GF(p) with p prime
+/// Returns order of x in GF(p) with p prime
 size_t order (int x, int p) {
     x = x % p;
     DAI_ASSERT( x != 0 );
 size_t order (int x, int p) {
     x = x % p;
     DAI_ASSERT( x != 0 );
@@ -289,7 +191,7 @@ size_t order (int x, int p) {
 }
 
 
 }
 
 
-// Returns whether n is a prime number
+/// Returns whether n is a prime number
 bool isPrime (size_t n) {
     bool result = true;
     for( size_t k = 2; (k < n) && result; k++ )
 bool isPrime (size_t n) {
     bool result = true;
     for( size_t k = 2; (k < n) && result; k++ )
@@ -299,8 +201,8 @@ bool isPrime (size_t n) {
 }
 
 
 }
 
 
-// Make a regular LDPC graph with N=6, j=2, K=4, k=3
-BipartiteGraph CreateSmallLDPCGraph() {
+/// Constructs 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;
 
     BipartiteGraph G;
     size_t N=4, j=3, K=4; // k=3;
 
@@ -319,13 +221,19 @@ BipartiteGraph CreateSmallLDPCGraph() {
 }
 
 
 }
 
 
-// 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 ) {
+/// Creates group-structured LDPC code
+/** 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)
+ *                      (p,j,k) = (37,3,4)
+ *                      (p,j,k) = (7,2,4)
+ *                      (p,j,k) = (29,2,4)
+ *
+ *  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;
     BipartiteGraph G;
 
     size_t n = j;
@@ -341,8 +249,8 @@ BipartiteGraph CreateGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) {
         if( order(b,p) == j )
             break;
     DAI_ASSERT( b != p );
         if( order(b,p) == j )
             break;
     DAI_ASSERT( b != p );
-//    cout << "# order(a=" << a << ") = " << order(a,p) << endl;
-//    cout << "# order(b=" << b << ") = " << order(b,p) << endl;
+    // cout << "# order(a=" << a << ") = " << order(a,p) << endl;
+    // cout << "# order(b=" << b << ") = " << order(b,p) << endl;
 
     DAI_ASSERT( N * n == K * k );
 
 
     DAI_ASSERT( N * n == K * k );
 
@@ -364,8 +272,8 @@ BipartiteGraph CreateGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) {
 }
 
 
 }
 
 
-// Make parity check table
-void MakeParityCheck( Real *result, size_t n, Real eps ) {
+// Constructs a parity check table
+void createParityCheck( Real *result, size_t n, Real eps ) {
     size_t N = 1 << n;
     for( size_t i = 0; i < N; i++ ) {
         size_t c = 0;
     size_t N = 1 << n;
     for( size_t i = 0; i < N; i++ ) {
         size_t c = 0;
@@ -381,116 +289,143 @@ void MakeParityCheck( Real *result, size_t n, Real eps ) {
 }
 
 
 }
 
 
-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";
+/// Predefined names of various factor graph types
+const char *FULL_TYPE        = "FULL";
+const char *DREG_TYPE        = "DREG";
+const char *LOOP_TYPE        = "LOOP";
+const char *TREE_TYPE        = "TREE";
+const char *GRID_TYPE        = "GRID";
+const char *GRID3D_TYPE      = "GRID3D";
+const char *HOI_TYPE         = "HOI";
+const char *LDPC_TYPE        = "LDPC";
+
+
+/// Possible LDPC structures
+DAI_ENUM(LDPCType,SMALL,GROUP,RANDOM);
 
 
 
 
+/// Main function
 int main( int argc, char *argv[] ) {
     try {
 int main( int argc, char *argv[] ) {
     try {
-        size_t N, K, k, d, j, n1, n2, n3;
-        size_t prime;
+        // Variables for storing command line arguments
         size_t seed;
         size_t seed;
-        Real beta, sigma_w, sigma_th, noise, mean_w, mean_th;
-        string type;
         size_t states = 2;
         size_t states = 2;
+        string type;
+        size_t d, N, K, k, j, n1, n2, n3, prime;
+        bool periodic = false;
+        FactorType ft;
+        LDPCType ldpc;
+        Real beta, sigma_w, sigma_th, mean_w, mean_th, noise;
 
         // Declare the supported options.
 
         // Declare the supported options.
-        po::options_description desc("Allowed options");
-        desc.add_options()
+        po::options_description opts("General command line options");
+        opts.add_options()
             ("help",     "produce help message")
             ("help",     "produce help message")
-            ("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')")
-            ("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<Real>(&beta),       "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)")
-            ("mean_w",   po::value<Real>(&mean_w),     "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
-            ("mean_th",  po::value<Real>(&mean_th),    "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
-            ("sigma_w",  po::value<Real>(&sigma_w),    "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
-            ("sigma_th", po::value<Real>(&sigma_th),   "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'")
-            ("noise",    po::value<Real>(&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')")
+            ("seed",     po::value<size_t>(&seed),   "random number seed (tries to read from /dev/urandom if not specified)")
+            ("states",   po::value<size_t>(&states), "number of states of each variable (default=2 for binary variables)")
+        ;
+
+        // Graph structure options
+        po::options_description opts_graph("Options for specifying graph structure");
+        opts_graph.add_options()
+            ("type",     po::value<string>(&type),   "factor graph type (one of 'FULL', 'DREG', 'LOOP', 'TREE', 'GRID', 'GRID3D', 'HOI', 'LDPC')")
+            ("d",        po::value<size_t>(&d),      "variable connectivity (only for type=='DREG');\n\t<d><N> should be even")
+            ("N",        po::value<size_t>(&N),      "number of variables (not for type=='GRID','GRID3D')")
+            ("n1",       po::value<size_t>(&n1),     "width of grid (only for type=='GRID','GRID3D')")
+            ("n2",       po::value<size_t>(&n2),     "height of grid (only for type=='GRID','GRID3D')")
+            ("n3",       po::value<size_t>(&n3),     "length of grid (only for type=='GRID3D')")
+            ("periodic", po::value<bool>(&periodic), "periodic grid? (only for type=='GRID','GRID3D'; default=0)")
+            ("K",        po::value<size_t>(&K),      "number of factors (only for type=='HOI','LDPC')")
+            ("k",        po::value<size_t>(&k),      "number of variables per factor (only for type=='HOI','LDPC')")
+        ;
+
+        // Factor options
+        po::options_description opts_factors("Options for specifying factors");
+        opts_factors.add_options()
+            ("factors",  po::value<FactorType>(&ft), "factor type (one of 'EXPGAUSS','POTTS','ISING')")
+            ("beta",     po::value<Real>(&beta),     "inverse temperature (ignored for factors=='ISING')")
+            ("mean_w",   po::value<Real>(&mean_w),   "mean of pairwise interactions w_{ij} (only for factors=='ISING')")
+            ("mean_th",  po::value<Real>(&mean_th),  "mean of unary interactions th_i (only for factors=='ISING')")
+            ("sigma_w",  po::value<Real>(&sigma_w),  "stddev of pairwise interactions w_{ij} (only for factors=='ISING')")
+            ("sigma_th", po::value<Real>(&sigma_th), "stddev of unary interactions th_i (only for factors=='ISING'")
+        ;
+
+        // LDPC options
+        po::options_description opts_ldpc("Options for specifying LDPC code factor graphs");
+        opts_ldpc.add_options()
+            ("ldpc",     po::value<LDPCType>(&ldpc), "type of LDPC code (one of 'SMALL','GROUP','RANDOM')")
+            ("j",        po::value<size_t>(&j),      "number of parity checks per bit (only for type=='LDPC')")
+            ("noise",    po::value<Real>(&noise),    "bitflip probability for binary symmetric channel (only for type=='LDPC')")
+            ("prime",    po::value<size_t>(&prime),  "prime number for construction of LDPC code (only for type=='LDPC' with ldpc='GROUP'))")
         ;
 
         ;
 
+        // All options
+        opts.add(opts_graph).add(opts_factors).add(opts_ldpc);
+
+        // Parse command line arguments
         po::variables_map vm;
         po::variables_map vm;
-        po::store(po::parse_command_line(argc, argv, desc), vm);
+        po::store(po::parse_command_line(argc, argv, opts), vm);
         po::notify(vm);
 
         po::notify(vm);
 
+        // Display help message if necessary
         if( vm.count("help") || !vm.count("type") ) {
         if( vm.count("help") || !vm.count("type") ) {
-            if( vm.count("type") ) {
-                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 << "<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 == 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', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d')" << 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;
+            cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
+            cout << "Usage: ./createfg [options]" << endl << endl;
+            cout << "Creates a factor graph according to the specified options." << endl << endl;
+            
+            cout << endl << opts << endl;
+            
+            cout << "The following factor graph types with pairwise interactions can be created:" << endl;
+            cout << "\t'FULL':   fully connected graph of <N> variables" << endl;
+            cout << "\t'DREG':   random regular graph of <N> variables where each variable is connected with <d> others" << endl;
+            cout << "\t'LOOP':   a single loop of <N> variables" << endl;
+            cout << "\t'TREE':   random tree-structured (acyclic, connected) graph of <N> variables" << endl;
+            cout << "\t'GRID':   2D grid of <n1>x<n2> variables" << endl;
+            cout << "\t'GRID3D': 3D grid of <n1>x<n2>x<n3> variables" << endl;
+            cout << "The following higher-order interactions factor graphs can be created:" << endl;
+            cout << "\t'HOI':    random factor graph consisting of <N> variables and <K> factors," << endl;
+            cout << "\t          each factor being an interaction of <k> variables." << endl;
+            cout << "The following LDPC code factor graphs can be created:" << endl;
+            cout << "\t'LDPC':   simulates LDPC decoding problem, using an LDPC code of <N> bits and <K>" << endl;
+            cout << "\t          parity checks, with <k> bits per check and <j> checks per bit, transmitted" << endl;
+            cout << "\t          on a binary symmetric channel with probability <noise> of flipping a bit." << endl;
+            cout << "\t          The transmitted codeword has all bits set to zero. The argument 'ldpc'" << endl;
+            cout << "\t          determines how the LDPC code is constructed: either using a group structure," << endl;
+            cout << "\t          or randomly, or a fixed small code with (N,K,k,j) = (4,4,3,3)." << endl << endl;
+
+            cout << "For all types except type=='LDPC', the factors have to be specified as well." << endl << endl;
+
+            cout << "EXPGAUSS factors (the default) are created by drawing all log-factor entries" << endl;
+            cout << "independently from a Gaussian with mean 0 and standard deviation <beta>." << endl << endl;
+            
+            cout << "In case of pairwise interactions, one can also choose POTTS factors, for which" << endl;
+            cout << "the log-factors are simply delta functions multiplied by the strength <beta>." << endl << endl;
+
+            cout << "For pairwise interactions and binary variables, one can also use ISING factors." << endl;
+            cout << "Here variables x1...xN are assumed to be +1/-1--valued, and unary interactions" << endl;
+            cout << "are of the form exp(th*xi) with th drawn from a Gaussian distribution with mean" << endl;
+            cout << "<mean_th> and standard deviation <sigma_th>, and pairwise interactions are of the" << endl;
+            cout << "form exp(w*xi*xj) with w drawn from a Gaussian distribution with mean <mean_w>" << endl;
+            cout << "and standard deviation <sigma_w>." << endl;
             return 1;
         }
 
             return 1;
         }
 
+        // Set default number of states
         if( !vm.count("states") )
             states = 2;
 
         if( !vm.count("states") )
             states = 2;
 
+        // Set default factor type
+        if( !vm.count("factors") )
+            ft = FactorType::EXPGAUSS;
+        // Check validness of factor type
+        if( ft == FactorType::POTTS )
+            if( type == HOI_TYPE )
+                throw "For factors=='POTTS', interactions should be pairwise (type!='HOI')";
+        if( ft == FactorType::ISING )
+            if( ((states != 2) || (type == HOI_TYPE)) )
+                throw "For factors=='ISING', variables should be binary (states==2) and interactions should be pairwise (type!='HOI')";
+
+        // Read random seed
         if( !vm.count("seed") ) {
             ifstream infile;
             bool success;
         if( !vm.count("seed") ) {
             ifstream infile;
             bool success;
@@ -506,128 +441,126 @@ int main( int argc, char *argv[] ) {
         }
         rnd_seed( seed );
 
         }
         rnd_seed( seed );
 
-        FactorGraph fg;
-
+        // Set default periodicity
+        if( !vm.count("periodic") )
+            periodic = false;
+
+        // Store some options in a PropertySet object
+        PropertySet options;
+        if( vm.count("mean_th") )
+            options.Set("mean_th", mean_th);
+        if( vm.count("sigma_th") )
+            options.Set("sigma_th", sigma_th);
+        if( vm.count("mean_w") )
+            options.Set("mean_w", mean_w);
+        if( vm.count("sigma_w") )
+            options.Set("sigma_w", sigma_w);
+        if( vm.count("beta") )
+            options.Set("beta", beta);
+
+        // Output some comments
         cout << "# Factor graph made by " << argv[0] << endl;
         cout << "# type = " << type << endl;
         cout << "# Factor graph made by " << argv[0] << endl;
         cout << "# type = " << type << endl;
+        cout << "# states = " << states << endl;
 
 
-        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, mean_w, mean_th, sigma_w, sigma_th, fg );
+        // The factor graph to be constructed
+        FactorGraph 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;
-        } else if( type == GRID_TYPE || type == GRID_TORUS_TYPE ) {
 #define NEED_ARG(name, desc) do { if(!vm.count(name)) throw "Please specify " desc " with --" name; } while(0);
 #define NEED_ARG(name, desc) do { if(!vm.count(name)) throw "Please specify " desc " with --" name; } while(0);
-            if( states > 2 ) {
-                NEED_ARG("N", "number of nodes");
+
+        if( type == FULL_TYPE || type == DREG_TYPE || type == LOOP_TYPE || type == TREE_TYPE || type == GRID_TYPE || type == GRID3D_TYPE ) {
+            // Pairwise interactions
+
+            // Check command line options
+            if( type == GRID_TYPE ) {
+                NEED_ARG("n1", "width of grid");
+                NEED_ARG("n2", "height of grid");
+                N = n1 * n2;
+            } else if( type == GRID3D_TYPE ) {
+                NEED_ARG("n1", "width of grid");
+                NEED_ARG("n2", "height of grid");
+                NEED_ARG("n3", "depth of grid");
+                N = n1 * n2 * n3;
+            } else
+                NEED_ARG("N", "number of variables");
+
+            if( states > 2 || ft == FactorType::POTTS ) {
                 NEED_ARG("beta", "stddev of log-factor entries");
             } else {
                 NEED_ARG("beta", "stddev of log-factor entries");
             } else {
-                NEED_ARG("N", "number of nodes");
                 NEED_ARG("mean_w", "mean of pairwise interactions");
                 NEED_ARG("mean_w", "mean of pairwise interactions");
-                NEED_ARG("mean_th", "mean of singleton interactions");
+                NEED_ARG("mean_th", "mean of unary interactions");
                 NEED_ARG("sigma_w", "stddev of pairwise interactions");
                 NEED_ARG("sigma_w", "stddev of pairwise interactions");
-                NEED_ARG("sigma_th", "stddev of singleton interactions");
+                NEED_ARG("sigma_th", "stddev of unary interactions");
             }
 
             }
 
-            size_t n = (size_t)sqrt((long double)N);
-            N = n * n;
-
-            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;
-
-            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;
+            if( type == DREG_TYPE )
+                NEED_ARG("d", "connectivity (number of neighboring variables of each variable)");
+
+            // Build pairwise interaction graph
+            GraphAL G;
+            if( type == FULL_TYPE )
+                G = createGraphFull( N );
+            else if( type == DREG_TYPE )
+                G = createGraphRegular( N, d );
+            else if( type == LOOP_TYPE )
+                G = createGraphLoop( N );
+            else if( type == TREE_TYPE )
+                G = createGraphTree( N );
+            else if( type == GRID_TYPE )
+                G = createGraphGrid( n1, n2, periodic );
+            else if( type == GRID3D_TYPE )
+                G = createGraphGrid3D( n1, n2, n3, periodic );
+
+            // Construct factor graph from pairwise interaction graph
+            fg = createFG( G, ft, states, options );
+
+            // Output some additional comments
+            if( type == GRID_TYPE || type == GRID3D_TYPE ) {
+                cout << "# n1 = " << n1 << endl;
+                cout << "# n2 = " << n2 << endl;
+                if( type == GRID3D_TYPE )
+                    cout << "# n3 = " << n3 << 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, 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( type == DREG_TYPE )
+                cout << "# d = " << d << endl;
+            cout << "# options = " << options << endl;
+        } else if( type == HOI_TYPE ) {
+            // Higher order interactions
 
 
-            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 );
+            // Check command line arguments
+            NEED_ARG("N", "number of variables");
+            NEED_ARG("K", "number of factors");
+            NEED_ARG("k", "number of variables per factor");
+            NEED_ARG("beta", "stddev of log-factor entries");
 
 
-            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;
-        } else if( type == HOI_TYPE ) {
-            if( !vm.count("N") || !vm.count("K") || !vm.count("k") || !vm.count("beta") )
-                throw "Please specify all required arguments";
+            // Create higher-order interactions factor graph
             do {
             do {
-                MakeHOIFG( N, K, k, beta, fg );
+                fg = createHOIFG( N, K, k, beta );
             } while( !fg.isConnected() );
 
             } while( !fg.isConnected() );
 
-            cout << "# N = " << N << endl;
+            // Output some additional comments
             cout << "# K = " << K << endl;
             cout << "# k = " << k << endl;
             cout << "# beta = " << beta << 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";
-
+        } else if( type == LDPC_TYPE ) {
+            // LDPC codes
+
+            // Check command line arguments
+            NEED_ARG("ldpc", "type of LDPC code");
+            NEED_ARG("noise", "bitflip probability for binary symmetric channel");
+
+            // Check more command line arguments (seperately for each LDPC type)
+            if( ldpc == LDPCType::RANDOM ) {
+                NEED_ARG("N", "number of variables");
+                NEED_ARG("K", "number of factors");
+                NEED_ARG("k", "number of variables per factor");
+                NEED_ARG("j", "number of parity checks per bit");
                 if( N * j != K * k )
                     throw "Parameters should satisfy N * j == K * k";
                 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";
+            } else if( ldpc == LDPCType::GROUP ) {
+                NEED_ARG("prime", "prime number");
+                NEED_ARG("k", "number of variables per factor");
+                NEED_ARG("j", "number of parity checks per bit");
 
                 if( !isPrime(prime) )
                     throw "Parameter <prime> should be prime";
 
                 if( !isPrime(prime) )
                     throw "Parameter <prime> should be prime";
@@ -638,41 +571,37 @@ int main( int argc, char *argv[] ) {
 
                 N = prime * k;
                 K = prime * j;
 
                 N = prime * k;
                 K = prime * j;
-            } else if( type == LDPC_SMALL_TYPE ) {
+            } else if( ldpc == LDPCType::SMALL ) {
                 N = 4;
                 K = 4;
                 j = 3;
                 k = 3;
             }
 
                 N = 4;
                 K = 4;
                 j = 3;
                 k = 3;
             }
 
+            // Output some additional comments
             cout << "# N = " << N << endl;
             cout << "# K = " << K << endl;
             cout << "# j = " << j << endl;
             cout << "# k = " << k << endl;
             cout << "# N = " << N << endl;
             cout << "# K = " << K << endl;
             cout << "# j = " << j << endl;
             cout << "# k = " << k << endl;
-            if( type == LDPC_GROUP_TYPE )
+            if( ldpc == LDPCType::GROUP )
                 cout << "# prime = " << prime << endl;
             cout << "# noise = " << noise << endl;
 
                 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
             Real likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise};
             Real *paritycheck = new Real[1 << k];
             // Construct likelihood and paritycheck factors
             Real likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise};
             Real *paritycheck = new Real[1 << k];
-            MakeParityCheck(paritycheck, k, 0.0);
+            createParityCheck(paritycheck, k, 0.0);
 
             // Create LDPC structure
             BipartiteGraph ldpcG;
             bool regular;
             do {
 
             // 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();
+                if( ldpc == LDPCType::GROUP )
+                    ldpcG = createGroupStructuredLDPCGraph( prime, j, k );
+                else if( ldpc == LDPCType::RANDOM )
+                    ldpcG = createRandomBipartiteGraph( N, K, j, k );
+                else if( ldpc == LDPCType::SMALL )
+                    ldpcG = createSmallLDPCGraph();
 
                 regular = true;
                 for( size_t i = 0; i < N; i++ )
 
                 regular = true;
                 for( size_t i = 0; i < N; i++ )
@@ -718,24 +647,17 @@ int main( int argc, char *argv[] ) {
 
             // Construct Factor Graph
             fg = FactorGraph( factors );
 
             // 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 {
+        } else
             throw "Invalid type";
             throw "Invalid type";
-        }
 
 
+        // Output additional comments
+        cout << "# N = " << fg.nrVars() << endl;
         cout << "# seed = " << seed << endl;
         cout << "# seed = " << seed << endl;
+
+        // Output factor graph
         cout << fg;
     } catch( const char *e ) {
         cout << fg;
     } catch( const char *e ) {
+        /// Display error message
         cerr << "Error: " << e << endl;
         return 1;
     }
         cerr << "Error: " << e << endl;
         return 1;
     }