Now uses GMP big integers to represent linear states / total number of states
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 11 Jul 2011 10:21:40 +0000 (12:21 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 11 Jul 2011 10:21:40 +0000 (12:21 +0200)
16 files changed:
ChangeLog
Makefile.CYGWIN
Makefile.LINUX
Makefile.MACOSX
Makefile.WINDOWS
include/dai/clustergraph.h
include/dai/doc.h
include/dai/factor.h
include/dai/index.h
include/dai/jtree.h
include/dai/util.h
include/dai/varset.h
src/factorgraph.cpp
src/jtree.cpp
tests/testbbp.cpp
utils/fginfo.cpp

index 07647dd..3a809ce 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+* Now uses GMP big integer type to represent linear states and the number of
+  states of VarSets (in some cases). This fixes a bug Jerome Maye found in the 
+  State() class in index.h. Unfortunately, the supplied Makefile.WINDOWS now
+  no longer works out-of-the-box.
 * Changed license from GPL v2+ to FreeBSD (BSD 2-clause)
 * Fixed numerical issues in MF, FBP and TRWBP (discovered in sparse branch)
 * Jerome Maye found a bug in the State() class in index.h; implemented
index f8e6491..14a3c9c 100644 (file)
@@ -43,7 +43,7 @@ CCINC=-Iinclude -I/cygdrive/e/cygwin/boost_1_42_0
 
 # LINKER
 # Standard libraries to include
-LIBS=-ldai
+LIBS=-ldai -lgmpxx -lgmp
 # For linking with BOOST libraries
 BOOSTLIBS_PO=-lboost_program_options
 BOOSTLIBS_UTF=-lboost_unit_test_framework
index 148a016..f11ca17 100644 (file)
@@ -44,7 +44,7 @@ CCINC=-Iinclude
 
 # LINKER
 # Standard libraries to include
-LIBS=-ldai
+LIBS=-ldai -lgmpxx -lgmp
 # For linking with BOOST libraries
 BOOSTLIBS_PO=-lboost_program_options-mt
 BOOSTLIBS_UTF=-lboost_unit_test_framework-mt
index 2963b13..dc46716 100644 (file)
@@ -41,7 +41,7 @@ CCINC=-Iinclude -I/opt/local/include
 
 # LINKER
 # Standard libraries to include
-LIBS=-ldai
+LIBS=-ldai -lgmpxx -lgmp
 # For linking with BOOST libraries
 BOOSTLIBS_PO=-lboost_program_options
 BOOSTLIBS_UTF=-lboost_unit_test_framework
index ba8323e..32f9edd 100644 (file)
@@ -8,8 +8,8 @@
 # This template contains configurations for compiling libDAI with Visual C++
 # under Windows (and GNU Make)
 #
-# It has been tested with Windows XP, Visual Studio 2008, MatLab R2008b and
-# Boost 1.42.0
+# TODO: the latest libDAI version depends on GMP (Windows users can use MPIR)
+# This file has to be updated in order to link with the right GMP/MPIR libraries
 #
 # To use it, simply copy this file to 'Makefile.conf' and adapt 'Makefile.conf'
 # to your local setup
index 1e2fe56..22886d3 100644 (file)
@@ -234,7 +234,7 @@ namespace dai {
                     varindices.insert( i );
 
                 // Do variable elimination
-                long double totalStates = 0.0;
+                BigInt totalStates = 0;
                 while( !varindices.empty() ) {
                     size_t i = f( cl, varindices );
                     VarSet Di = cl.elimVar( i );
index 7ed4a11..24a9aab 100644 (file)
  *  \section compatibility Compatibility
  *  
  *  The code has been developed under Debian GNU/Linux with the GCC compiler suite.
- *  libDAI compiles successfully with g++ versions 3.4 up to 4.4.
+ *  libDAI compiles successfully with g++ versions 3.4 up to 4.6.
  *
  *  libDAI has also been successfully compiled with MS Visual Studio 2008 under Windows
  *  (but not all build targets are supported yet) and with Cygwin under Windows.
  *    - GNU make
  *    - recent boost C++ libraries (at least version 1.37; however,
  *      version 1.37 shipped with Ubuntu 9.04 is known not to work)
+ *    - GMP library (or the Windows port called MPIR)
  *    - doxygen (only for building the documentation)
  *    - graphviz (only for using some of the libDAI command line utilities)
  *    - CImg library (only for building the image segmentation example)
  * 
  *  On Debian/Ubuntu, you can easily install the required packages with a single command:
- *  <pre>  apt-get install g++ make doxygen graphviz libboost-dev libboost-graph-dev libboost-program-options-dev libboost-test-dev cimg-dev</pre>
+ *  <pre>  apt-get install g++ make doxygen graphviz libboost-dev libboost-graph-dev libboost-program-options-dev libboost-test-dev libgmp-dev cimg-dev</pre>
  *  (root permissions needed).
  *
  *  On Mac OS X (10.4 is known to work), these packages can be installed easily via MacPorts.
  *  If MacPorts is not already installed, install it according to the instructions at http://www.macports.org/.
  *  Then, a simple 
- *    <pre>  sudo port install gmake boost doxygen graphviz</pre>
+ *    <pre>  sudo port install gmake boost gmp doxygen graphviz</pre>
  *  should be enough to install everything that is needed.
  *  
  *  On Cygwin, the prebuilt Cygwin package boost-1.33.1-x is known not to work.
  *  You need:
  *  - A recent version of MicroSoft Visual Studio (2008 is known to work)
  *  - recent boost C++ libraries (version 1.37 or higher)
+ *  - GMP or MPIR library
  *  - GNU make (can be obtained from http://gnuwin32.sourceforge.net)
  *  - CImg library (only for building the image segmentation example)
  *
  *    <pre>
  *    bjam --with-graph --with-math --with-program_options --with-test link=static runtime-link=shared</pre>
  *
+ *  \subsection build-windows-gmp Building GMP or MPIR under Windows
+ *  
+ *  Information about how to build GPR or MPIR under Windows can be found on the internet.
+ *  The user has to update Makefile.WINDOWS in order to link with the GPR/MPIR libraries.
+ *
  *  \subsection build-windows-libdai Building libDAI
  *
  *  To build the source, copy Makefile.WINDOWS to Makefile.conf. Then, edit 
index 01079be..f4d713e 100644 (file)
@@ -69,13 +69,13 @@ class TFactor {
         TFactor( const Var &v ) : _vs(v), _p(v.states()) {}
 
         /// Constructs factor depending on variables in \a vars with uniform distribution
-        TFactor( const VarSet& vars ) : _vs(vars), _p((size_t)_vs.nrStates()) {
-            DAI_ASSERT( _vs.nrStates() <= std::numeric_limits<std::size_t>::max() );
+        TFactor( const VarSet& vars ) : _vs(vars), _p() {
+            _p = TProb<T>( BigInt_size_t( _vs.nrStates() ) );
         }
 
         /// Constructs factor depending on variables in \a vars with all values set to \a p
-        TFactor( const VarSet& vars, T p ) : _vs(vars), _p((size_t)_vs.nrStates(),p) {
-            DAI_ASSERT( _vs.nrStates() <= std::numeric_limits<std::size_t>::max() );
+        TFactor( const VarSet& vars, T p ) : _vs(vars), _p() {
+            _p = TProb<T>( BigInt_size_t( _vs.nrStates() ), p );
         }
 
         /// Constructs factor depending on variables in \a vars, copying the values from a std::vector<>
@@ -93,8 +93,9 @@ class TFactor {
         /** \param vars contains the variables that the new factor should depend on.
          *  \param p Points to array of values to be added.
          */
-        TFactor( const VarSet& vars, const T* p ) : _vs(vars), _p(p, p + (size_t)_vs.nrStates(), (size_t)_vs.nrStates()) {
-            DAI_ASSERT( _vs.nrStates() <= std::numeric_limits<std::size_t>::max() );
+        TFactor( const VarSet& vars, const T* p ) : _vs(vars), _p() {
+            size_t N = BigInt_size_t( _vs.nrStates() );
+            _p = TProb<T>( p, p + N, N );
         }
 
         /// Constructs factor depending on variables in \a vars, copying the values from \a p
@@ -104,7 +105,7 @@ class TFactor {
 
         /// Constructs factor depending on variables in \a vars, permuting the values given in \a p accordingly
         TFactor( const std::vector<Var> &vars, const std::vector<T> &p ) : _vs(vars.begin(), vars.end(), vars.size()), _p(p.size()) {
-            size_t nrStates = 1;
+            BigInt nrStates = 1;
             for( size_t i = 0; i < vars.size(); i++ )
                 nrStates *= vars[i].states();
             DAI_ASSERT( nrStates == p.size() );
@@ -344,8 +345,7 @@ class TFactor {
             else {
                 TFactor<T> f(*this); // make a copy
                 _vs |= g._vs;
-                DAI_ASSERT( _vs.nrStates() < std::numeric_limits<std::size_t>::max() );
-                size_t N = (size_t)_vs.nrStates();
+                size_t N = BigInt_size_t( _vs.nrStates() );
 
                 IndexFor i_f( f._vs, _vs );
                 IndexFor i_g( g._vs, _vs );
@@ -403,8 +403,7 @@ class TFactor {
                 result._p = _p.pwBinaryTr( g._p, op );
             } else {
                 result._vs = _vs | g._vs;
-                DAI_ASSERT( result._vs.nrStates() < std::numeric_limits<std::size_t>::max() );
-                size_t N = (size_t)result._vs.nrStates();
+                size_t N = BigInt_size_t( result._vs.nrStates() );
 
                 IndexFor i_f( _vs, result.vars() );
                 IndexFor i_g( g._vs, result.vars() );
index afff672..911a03d 100644 (file)
@@ -346,15 +346,8 @@ class multifor {
  *
  *  \see dai::calcLinearState(), dai::calcState()
  *
- *  \idea Make the State class a more prominent part of libDAI 
- *  (and document it clearly, explaining the concept of state); 
- *  add more optimized variants of the State class like IndexFor 
- *  (e.g. for TFactor<>::slice()).
- *
- *  \todo The State class is dangerous because currently it represents a state both as
- *  an integer and as a map<Var, size_t>; if the cardinality of the joint state space
- *  becomes too large, the state no longer fits into an integer, and the two representations
- *  are no longer consistent... this should be fixed!
+ *  \idea Make the State class a more prominent part of libDAI (and document it clearly, explaining the concept of state); 
+ *  add more optimized variants of the State class like IndexFor (e.g. for TFactor<>::slice()).
  */
 class State {
     private:
@@ -362,7 +355,7 @@ class State {
         typedef std::map<Var, size_t> states_type;
 
         /// Current state (represented linearly)
-        long                          state;
+        BigInt                        state;
 
         /// Current state (represented as a map)
         states_type                   states;
@@ -372,13 +365,13 @@ class State {
         State() : state(0), states() {}
 
         /// Construct from VarSet \a vs and corresponding linear state \a linearState
-        State( const VarSet &vs, size_t linearState=0 ) : state(linearState), states() {
+        State( const VarSet &vs, BigInt linearState=0 ) : state(linearState), states() {
             if( linearState == 0 )
                 for( VarSet::const_iterator v = vs.begin(); v != vs.end(); v++ )
                     states[*v] = 0;
             else {
                 for( VarSet::const_iterator v = vs.begin(); v != vs.end(); v++ ) {
-                    states[*v] = linearState % v->states();
+                    states[*v] = BigInt_size_t( linearState % v->states() );
                     linearState /= v->states();
                 }
                 DAI_ASSERT( linearState == 0 );
@@ -402,7 +395,7 @@ class State {
         /// Return current linear state
         operator size_t() const {
             DAI_ASSERT( valid() );
-            return( state );
+            return( BigInt_size_t( state ) );
         }
 
         /// Inserts a range of variable-state pairs, changing the current state
@@ -432,9 +425,9 @@ class State {
         }
 
         /// Return linear state of variables in \a vs, assuming that variables that are not in \c *this are in state 0
-        size_t operator() ( const VarSet &vs ) const {
-            size_t vs_state = 0;
-            size_t prod = 1;
+        BigInt operator() ( const VarSet &vs ) const {
+            BigInt vs_state = 0;
+            BigInt prod = 1;
             for( VarSet::const_iterator v = vs.begin(); v != vs.end(); v++ ) {
                 states_type::const_iterator entry = states.find( *v );
                 if( entry != states.end() )
index 1bc98e2..d52d05b 100644 (file)
@@ -201,7 +201,7 @@ class JTree : public DAIAlgRG {
  *  \throws OUT_OF_MEMORY if the total number of states becomes larger than maxStates
  *  \return a pair (number of variables in largest clique, number of states in largest clique)
  */
-std::pair<size_t,long double> boundTreewidth( const FactorGraph &fg, greedyVariableElimination::eliminationCostFunction fn, size_t maxStates=0 );
+std::pair<size_t,BigInt> boundTreewidth( const FactorGraph &fg, greedyVariableElimination::eliminationCostFunction fn, size_t maxStates=0 );
 
 
 } // end of namespace dai
index e8ecdfc..8ca3073 100644 (file)
@@ -24,6 +24,7 @@
 #include <boost/lexical_cast.hpp>
 #include <algorithm>
 #include <cerrno>
+#include <gmpxx.h>
 
 #include <dai/exceptions.h>
 
@@ -79,6 +80,15 @@ namespace dai {
 /// Real number (alias for \c double, which could be changed to <tt>long double</tt> if necessary)
 typedef double Real;
 
+/// Arbitrary precision integer number
+typedef mpz_class BigInt;
+
+/// Safe down-cast of big integer to size_t
+inline size_t BigInt_size_t( const BigInt &N ) {
+    DAI_ASSERT( N <= std::numeric_limits<std::size_t>::max() );
+    return N.get_ui();
+}
+
 /// Returns true if argument is NAN (Not A Number)
 bool isnan( Real x );
 
index d82ea5a..d1a76ae 100644 (file)
@@ -127,8 +127,8 @@ class VarSet : public SmallSet<Var> {
          *  number of possible values ("states") of variable \f$x_l\f$, the number of
          *  joint configurations of the variables in \f$\{x_l\}_{l\in L}\f$ is given by \f$\prod_{l\in L} S_l\f$.
          */
-        long double nrStates() const {
-            long double states = 1.0;
+        BigInt nrStates() const {
+            BigInt states = 1;
             for( VarSet::const_iterator n = begin(); n != end(); n++ )
                 states *= n->states();
             return states;
index 7581713..2094931 100644 (file)
@@ -348,7 +348,7 @@ Real FactorGraph::logScore( const std::vector<size_t>& statevec ) const {
     // by summing the log factor entries of the factors that correspond to this joint configuration
     Real lS = 0.0;
     for( size_t I = 0; I < nrFactors(); I++ )
-        lS += dai::log( factor(I)[S(factor(I).vars())] );
+        lS += dai::log( factor(I)[BigInt_size_t(S(factor(I).vars()))] );
     return lS;
 }
 
index b517bb3..4fb5e42 100644 (file)
@@ -96,7 +96,7 @@ JTree::JTree( const FactorGraph &fg, const PropertySet &opts, bool automatic ) :
             cerr << "VarElim result: " << ElimVec << endl;
 
         // Estimate memory needed (rough upper bound)
-        long double memneeded = 0;
+        BigInt memneeded = 0;
         foreach( const VarSet& cl, ElimVec )
             memneeded += cl.nrStates();
         memneeded *= sizeof(Real) * fudge;
@@ -393,10 +393,10 @@ Real JTree::logZ() const {
 
 size_t JTree::findEfficientTree( const VarSet& vs, RootedTree &Tree, size_t PreviousRoot ) const {
     // find new root clique (the one with maximal statespace overlap with vs)
-    long double maxval = 0.0;
+    BigInt maxval = 0;
     size_t maxalpha = 0;
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        long double val = VarSet(vs & OR(alpha).vars()).nrStates();
+        BigInt val = VarSet(vs & OR(alpha).vars()).nrStates();
         if( val > maxval ) {
             maxval = val;
             maxalpha = alpha;
@@ -556,7 +556,7 @@ Factor JTree::calcMarginal( const VarSet& vs ) {
 }
 
 
-std::pair<size_t,long double> boundTreewidth( const FactorGraph &fg, greedyVariableElimination::eliminationCostFunction fn, size_t maxStates ) {
+std::pair<size_t,BigInt> boundTreewidth( const FactorGraph &fg, greedyVariableElimination::eliminationCostFunction fn, size_t maxStates ) {
     // Create cluster graph from factor graph
     ClusterGraph _cg( fg, true );
 
@@ -565,11 +565,11 @@ std::pair<size_t,long double> boundTreewidth( const FactorGraph &fg, greedyVaria
 
     // Calculate treewidth
     size_t treewidth = 0;
-    double nrstates = 0.0;
+    BigInt nrstates = 0.0;
     for( size_t i = 0; i < ElimVec.size(); i++ ) {
         if( ElimVec[i].size() > treewidth )
             treewidth = ElimVec[i].size();
-        long double s = ElimVec[i].nrStates();
+        BigInt s = ElimVec[i].nrStates();
         if( s > nrstates )
             nrstates = s;
     }
index a06f5ff..168e5d1 100644 (file)
@@ -31,7 +31,6 @@ int main( int argc, char *argv[] ) {
         Real   tol = 1.0e-9;
         size_t maxiter = 10000;
         Real   damping = 0.0;
-        BBP::Properties::UpdateType updates = BBP::Properties::UpdateType::PAR;
 
         // Store the constants in a PropertySet object
         PropertySet opts;
index efdd899..3318179 100644 (file)
@@ -95,7 +95,7 @@ int main( int argc, char *argv[] ) {
         cout << "Pairwise interactions? " << fg.isPairwise() << endl;
         
         // Calculate treewidth using various heuristics
-        std::pair<size_t,double> tw;
+        std::pair<size_t,BigInt> tw;
         cout << "Treewidth (MinNeighbors):     ";
         try {
             tw = boundTreewidth(fg, &eliminationCost_MinNeighbors, maxstates );
@@ -141,7 +141,7 @@ int main( int argc, char *argv[] ) {
         }
         
         // Calculate total state space
-        long double stsp = 1.0;
+        BigInt stsp = 1;
         for( size_t i = 0; i < fg.nrVars(); i++ )
             stsp *= fg.var(i).states();
         cout << "Total state space:   " << stsp << endl;
@@ -149,9 +149,9 @@ int main( int argc, char *argv[] ) {
         cout << "Type: " << (fg.isPairwise() ? "pairwise" : "higher order") << " interactions, " << (fg.isBinary() ? "binary" : "nonbinary") << " variables" << endl;
 
         // Calculate complexity for LCBP
-        long double cavsum_lcbp = 0.0;
-        long double cavsum_lcbp2 = 0.0;
-        long double max_Delta_size = 0.0;
+        BigInt cavsum_lcbp = 0;
+        BigInt cavsum_lcbp2 = 0;
+        BigInt max_Delta_size = 0;
         map<size_t,size_t> cavsizes;
         for( size_t i = 0; i < fg.nrVars(); i++ ) {
             VarSet di = fg.delta(i);
@@ -159,7 +159,7 @@ int main( int argc, char *argv[] ) {
                 cavsizes[di.size()]++;
             else
                 cavsizes[di.size()] = 1;
-            long double Ds = fg.Delta(i).nrStates();
+            BigInt Ds = fg.Delta(i).nrStates();
             if( Ds > max_Delta_size )
                 max_Delta_size = Ds;
             cavsum_lcbp += di.nrStates();