Merge branch 'eaton'

Conflicts:

	Makefile
	src/bp.cpp
	tests/testdai.cpp

Summary of changes done in branch 'eaton':
* [Frederik Eaton] Improved parsing code in tests/testdai to allow recursive
  alias substitutions
* Interface changes:
    TProb<T>::minVal()     -> TProb<T>::min()
    TFactor<T>::minVal()   -> TFactor<T>::min()
    TProb<T>::maxVal()     -> TProb<T>::max()
    TFactor<T>::maxVal()   -> TFactor<T>::max()
    TProb<T>::totalSum()   -> TProb<T>sum()
    TFactor<T>::totalSum() -> TFactor<T>::sum()
* [Frederik Eaton] Added methods setUniform (), sumAbs(), argmax() to TProb<T>
* [Frederik Eaton] Added TAGS, lib targets to Makefile
* [Frederik Eaton] Add useful aliases to aliases.conf

diff --git a/ChangeLog b/ChangeLog
index 9b1c5a5..c1f8a37 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,15 @@
+* [Frederik Eaton] Improved parsing code in tests/testdai to allow recursive
+  alias substitutions
+* Interface changes:
+    TProb<T>::minVal()     -> TProb<T>::min()
+    TFactor<T>::minVal()   -> TFactor<T>::min()
+    TProb<T>::maxVal()     -> TProb<T>::max()
+    TFactor<T>::maxVal()   -> TFactor<T>::max()
+    TProb<T>::totalSum()   -> TProb<T>sum()
+    TFactor<T>::totalSum() -> TFactor<T>::sum()
+* [Frederik Eaton] Added methods setUniform (), sumAbs(), argmax() to TProb<T>
+* [Frederik Eaton] Added TAGS, lib targets to Makefile
+* [Frederik Eaton] Add useful aliases to aliases.conf
 * [Frederik Eaton] Removed unnecessary copy constructors and assignment
   operators
 * [Frederik Eaton] Change cout to cerr in warnings and error messages

diff --git a/Makefile b/Makefile
index 0c71239..2d63911 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -85,7 +85,7 @@ ifeq ($(OS),MACOSX)
   BOOSTLIBS=-lboost_program_options-mt
 endif
 ifeq ($(OS),LINUX)
-  BOOSTLIBS=-lboost_program_options
+  BOOSTLIBS=-lboost_program_options-mt
 endif
 ifeq ($(OS),WINDOWS)
   BOOSTLIBS=/LIBPATH:C:\boost_1_36_0\stage\lib
@@ -125,7 +125,7 @@ ifeq ($(OS),MACOSX)
 endif
 
 # Build targets
-TARGETS=tests utils $(LIB)/libdai$(LE) examples testregression
+TARGETS=tests utils lib examples testregression
 ifneq ($(OS),WINDOWS)
   TARGETS:=$(TARGETS) doc 
 endif
@@ -216,6 +216,8 @@ tests : tests/testdai$(EE)
 
 utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE)
 
+lib: $(LIB)/libdai$(LE) 
+
 
 # OBJECTS
 ##########
@@ -367,6 +369,10 @@ endif
 doc : $(INC)/*.h $(SRC)/*.cpp examples/*.cpp doxygen.conf
        doxygen doxygen.conf
 
+TAGS:
+       etags src/*.cpp include/dai/*.h tests/*.cpp utils/*.cpp
+       ctags src/*.cpp include/dai/*.h tests/*.cpp utils/*.cpp
+
 
 # CLEAN
 ########

diff --git a/include/dai/factor.h b/include/dai/factor.h
index 4f91f14..b441116 100644 (file)
--- a/include/dai/factor.h
+++ b/include/dai/factor.h
@@ -393,16 +393,16 @@ template <typename T> class TFactor {
         bool hasNegatives() const { return _p.hasNegatives(); }
 
         /// Returns total sum of values
-        T totalSum() const { return _p.totalSum(); }
+        T sum() const { return _p.sum(); }
 
         /// Returns maximum absolute value
         T maxAbs() const { return _p.maxAbs(); }
 
         /// Returns maximum value
-        T maxVal() const { return _p.maxVal(); }
+        T max() const { return _p.max(); }
 
         /// Returns minimum value
-        T minVal() const { return _p.minVal(); }
+        T min() const { return _p.min(); }
 
         /// Returns entropy of *this TFactor
         Real entropy() const { return _p.entropy(); }
@@ -486,8 +486,8 @@ template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const
                                 bs = i.states();
                             else
                                 as = j.states();
-                            T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).maxVal();
-                            T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).maxVal();
+                            T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).max();
+                            T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).max();
                             T f = f1 * f2;
                             if( f > max )
                                 max = f;

diff --git a/include/dai/prob.h b/include/dai/prob.h
index 9a8148d..238f32f 100644 (file)
--- a/include/dai/prob.h
+++ b/include/dai/prob.h
@@ -150,6 +150,12 @@ template <typename T> class TProb {
                     _p[i] = 0;
             return *this;
         }
+        
+        /// Set all entries to 1.0/size()
+        TProb<T>& setUniform () {
+            fill(1.0/size());
+            return *this;
+        }
 
         /// Sets entries that are smaller than epsilon to epsilon
         TProb<T>& makePositive( Real epsilon ) {
@@ -403,11 +409,19 @@ template <typename T> class TProb {
         }
 
         /// Returns sum of all entries
-        T totalSum() const {
+        T sum() const {
             T Z = std::accumulate( _p.begin(),  _p.end(), (T)0 );
             return Z;
         }
 
+        /// Return sum of absolute value of all entries
+        T sumAbs() const {
+            T s = 0;
+            for( size_t i = 0; i < size(); i++ )
+                s += fabs( (Real) _p[i] );
+            return s;
+        }
+
         /// Returns maximum absolute value of all entries
         T maxAbs() const {
             T Z = 0;
@@ -420,22 +434,35 @@ template <typename T> class TProb {
         }
 
         /// Returns maximum value of all entries
-        T maxVal() const {
+        T max() const {
             T Z = *std::max_element( _p.begin(), _p.end() );
             return Z;
         }
 
         /// Returns minimum value of all entries
-        T minVal() const {
+        T min() const {
             T Z = *std::min_element( _p.begin(), _p.end() );
             return Z;
         }
 
+        /// Returns {arg,}maximum value
+        std::pair<size_t,T> argmax() const {
+            T max = _p[0];
+            size_t arg = 0;
+            for( size_t i = 1; i < size(); i++ ) {
+              if( _p[i] > max ) {
+                max = _p[i];
+                arg = i;
+              }
+            }
+            return make_pair(arg,max);
+        }
+
         /// Normalizes vector using the specified norm
         T normalize( NormType norm=NORMPROB ) {
             T Z = 0.0;
             if( norm == NORMPROB )
-                Z = totalSum();
+                Z = sum();
             else if( norm == NORMLINF )
                 Z = maxAbs();
             if( Z == 0.0 )
@@ -478,7 +505,7 @@ template <typename T> class TProb {
 
         /// Returns a random index, according to the (normalized) distribution described by *this
         size_t draw() {
-            double x = rnd_uniform() * totalSum();
+            double x = rnd_uniform() * sum();
             T s = 0;
             for( size_t i = 0; i < size(); i++ ) {
                 s += _p[i];

diff --git a/include/dai/properties.h b/include/dai/properties.h
index 7b94a07..196fcd4 100644 (file)
--- a/include/dai/properties.h
+++ b/include/dai/properties.h
@@ -36,6 +36,7 @@
 #include <cassert>
 #include <typeinfo>
 #include <dai/exceptions.h>
+#include <dai/util.h>
 
 
 namespace dai {
@@ -72,6 +73,15 @@ class PropertySet : private std::map<PropertyKey, PropertyValue> {
         /// Sets a property
         PropertySet & Set(const PropertyKey &key, const PropertyValue &val) { this->operator[](key) = val; return *this; }
 
+        /// Set properties according to those in newProps, overriding properties that already exist with new values
+        PropertySet & Set( const PropertySet &newProps ) {
+            const std::map<PropertyKey, PropertyValue> *m = &newProps;
+            foreach(value_type i, *m) {
+                Set(i.first, i.second);
+            }
+            return *this;
+        }
+
         /// Gets a property, casted as ValueType
         template<typename ValueType>
         ValueType GetAs(const PropertyKey &key) const {

diff --git a/src/bp.cpp b/src/bp.cpp
index e83c65b..0b0046e 100644 (file)
--- a/src/bp.cpp
+++ b/src/bp.cpp
@@ -212,7 +212,7 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
             }
         }
         if( props.logdomain ) {
-            prod -= prod.maxVal();
+            prod -= prod.max();
             prod.takeExp();
         }
 
@@ -396,7 +396,7 @@ Factor BP::beliefV( size_t i ) const {
     calcBeliefV( i, p );
 
     if( props.logdomain ) {
-        p -= p.maxVal();
+        p -= p.max();
         p.takeExp();
     }
 
@@ -453,7 +453,7 @@ Factor BP::beliefF( size_t I ) const {
         calcBeliefF( I, prod );
 
         if( props.logdomain ) {
-            prod -= prod.maxVal();
+            prod -= prod.max();
             prod.takeExp();
         }
         prod.normalize();

diff --git a/src/exactinf.cpp b/src/exactinf.cpp
index 1b31aba..b4b970a 100644 (file)
--- a/src/exactinf.cpp
+++ b/src/exactinf.cpp
@@ -86,7 +86,7 @@ double ExactInf::run() {
     for( size_t I = 0; I < nrFactors(); I++ )
         P *= factor(I);
 
-    Real Z = P.totalSum();
+    Real Z = P.sum();
     _logZ = std::log(Z);
     for( size_t i = 0; i < nrVars(); i++ )
         _beliefsV[i] = P.marginal(var(i));

diff --git a/src/hak.cpp b/src/hak.cpp
index 40e34ee..250e9d9 100644 (file)
--- a/src/hak.cpp
+++ b/src/hak.cpp
@@ -473,14 +473,14 @@ vector<Factor> HAK::beliefs() const {
 
 
 Real HAK::logZ() const {
-    Real sum = 0.0;
+    Real s = 0.0;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum += IR(beta).c() * Qb(beta).entropy();
+        s += IR(beta).c() * Qb(beta).entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        sum += OR(alpha).c() * Qa(alpha).entropy();
-        sum += (OR(alpha).log(true) * Qa(alpha)).totalSum();
+        s += OR(alpha).c() * Qa(alpha).entropy();
+        s += (OR(alpha).log(true) * Qa(alpha)).sum();
     }
-    return sum;
+    return s;
 }
 
 

diff --git a/src/jtree.cpp b/src/jtree.cpp
index 177b405..d1bfceb 100644 (file)
--- a/src/jtree.cpp
+++ b/src/jtree.cpp
@@ -308,14 +308,14 @@ double JTree::run() {
 
 
 Real JTree::logZ() const {
-    Real sum = 0.0;
+    Real s = 0.0;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum += IR(beta).c() * Qb[beta].entropy();
+        s += IR(beta).c() * Qb[beta].entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        sum += OR(alpha).c() * Qa[alpha].entropy();
-        sum += (OR(alpha).log(true) * Qa[alpha]).totalSum();
+        s += OR(alpha).c() * Qa[alpha].entropy();
+        s += (OR(alpha).log(true) * Qa[alpha]).sum();
     }
-    return sum;
+    return s;
 }
 
 

diff --git a/src/mf.cpp b/src/mf.cpp
index baf7835..83b98e3 100644 (file)
--- a/src/mf.cpp
+++ b/src/mf.cpp
@@ -189,10 +189,10 @@ vector<Factor> MF::beliefs() const {
 
 
 Real MF::logZ() const {
-    Real sum = 0.0;
+    Real s = 0.0;
     
     for(size_t i=0; i < nrVars(); i++ )
-        sum -= beliefV(i).entropy();
+        s -= beliefV(i).entropy();
     for(size_t I=0; I < nrFactors(); I++ ) {
         Factor henk;
         foreach( const Neighbor &j, nbF(I) )  // for all j in I
@@ -201,10 +201,10 @@ Real MF::logZ() const {
         Factor piet;
         piet = factor(I).log(true);
         piet *= henk;
-        sum -= piet.totalSum();
+        s -= piet.sum();
     }
 
-    return -sum;
+    return -s;
 }
 
 

diff --git a/src/treeep.cpp b/src/treeep.cpp
index c9e45e3..2c88f98 100644 (file)
--- a/src/treeep.cpp
+++ b/src/treeep.cpp
@@ -177,23 +177,23 @@ void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<F
     // Normalize Qa and Qb
     _logZ = 0.0;
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ ) {
-        _logZ += log(Qa[_a[alpha]].totalSum());
+        _logZ += log(Qa[_a[alpha]].sum());
         Qa[_a[alpha]].normalize();
     }
     for( size_t beta = 0; beta < _Qb.size(); beta++ ) {
-        _logZ -= log(Qb[_b[beta]].totalSum());
+        _logZ -= log(Qb[_b[beta]].sum());
         Qb[_b[beta]].normalize();
     }
 }
 
 
 double TreeEP::TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const {
-    double sum = 0.0;
+    double s = 0.0;
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
-        sum += (Qa[_a[alpha]] * _Qa[alpha].log(true)).totalSum();
+        s += (Qa[_a[alpha]] * _Qa[alpha].log(true)).sum();
     for( size_t beta = 0; beta < _Qb.size(); beta++ )
-        sum -= (Qb[_b[beta]] * _Qb[beta].log(true)).totalSum();
-    return sum + _logZ;
+        s -= (Qb[_b[beta]] * _Qb[beta].log(true)).sum();
+    return s + _logZ;
 }
 
 
@@ -431,24 +431,24 @@ double TreeEP::run() {
 
 
 Real TreeEP::logZ() const {
-    double sum = 0.0;
+    double s = 0.0;
 
     // entropy of the tree
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum -= Qb[beta].entropy();
+        s -= Qb[beta].entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        sum += Qa[alpha].entropy();
+        s += Qa[alpha].entropy();
 
     // energy of the on-tree factors
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        sum += (OR(alpha).log(true) * Qa[alpha]).totalSum();
+        s += (OR(alpha).log(true) * Qa[alpha]).sum();
 
     // energy of the off-tree factors
     for( size_t I = 0; I < nrFactors(); I++ )
         if( offtree(I) )
-            sum += (_Q.find(I))->second.logZ( Qa, Qb );
+            s += (_Q.find(I))->second.logZ( Qa, Qb );
     
-    return sum;
+    return s;
 }
 
 

diff --git a/tests/aliases.conf b/tests/aliases.conf
index d25cb74..f031109 100644 (file)
--- a/tests/aliases.conf
+++ b/tests/aliases.conf
@@ -83,5 +83,20 @@ LCBP_PAIR2CAV_SEQRND:           LC[cavity=PAIR2,reinit=0,updates=SEQRND,maxiter=
 LCBP_PAIR2CAV_NONE:             LC[cavity=PAIR2,reinit=0,updates=SEQFIX,maxiter=0,cavainame=BP,cavaiopts=[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0],tol=1e-9,verbose=0]
 LCBP_UNICAV_SEQFIX:             LC[cavity=UNIFORM,updates=SEQFIX,maxiter=10000,tol=1e-9,verbose=0,cavaiopts=[],cavainame=NONE]
 LCBP_UNICAV_SEQRND:             LC[cavity=UNIFORM,updates=SEQRND,maxiter=10000,tol=1e-9,verbose=0,cavaiopts=[],cavainame=NONE]
+
 LCTREEEP:                       LC[cavity=FULL,reinit=1,updates=SEQFIX,maxiter=10000,cavainame=TREEEP,cavaiopts=[type=ORG,tol=1e-9,maxiter=10000,verbose=0],tol=1e-9,verbose=0]
 LCMF:                           LC[cavity=FULL,reinit=1,updates=SEQFIX,maxiter=10000,cavainame=MF,cavaiopts=[tol=1e-9,maxiter=10000,verbose=0],tol=1e-9,verbose=0]
+LCBP:                           LCBP_FULLCAVin_SEQRND
+
+# --- GIBBS -------------------
+
+GIBBS:                          GIBBS[iters=1000,verbose=0]
+GIBBS_1e1:                      GIBBS[iters=10]
+GIBBS_1e2:                      GIBBS[iters=100]
+GIBBS_1e3:                      GIBBS[iters=1000]
+GIBBS_1e4:                      GIBBS[iters=10000]
+GIBBS_1e5:                      GIBBS[iters=100000]
+GIBBS_1e6:                      GIBBS[iters=1000000]
+GIBBS_1e7:                      GIBBS[iters=10000000]
+GIBBS_1e8:                      GIBBS[iters=100000000]
+GIBBS_1e9:                      GIBBS[iters=1000000000]

diff --git a/tests/testdai.cpp b/tests/testdai.cpp
index 0774836..b79b5da 100644 (file)
--- a/tests/testdai.cpp
+++ b/tests/testdai.cpp
@@ -76,7 +76,7 @@ class TestDAI {
                 delete obj;
         }
 
-        string identify() { 
+        string identify() const {
             if( obj != NULL )
                 return obj->identify(); 
             else
@@ -155,47 +155,51 @@ class TestDAI {
 };
 
 
-pair<string, PropertySet> parseMethod( const string &_s, const map<string,string> & aliases ) {
-    // s = first part of _s, until '['
-    string::size_type pos = _s.find_first_of('[');
-    string s;
-    if( pos == string::npos )
-        s = _s;
-    else
-        s = _s.substr(0,pos);
-
-    // if the first part is an alias, substitute
-    if( aliases.find(s) != aliases.end() )
-        s = aliases.find(s)->second;
+pair<string, PropertySet> parseMethodRaw( const string &s ) {
+    string::size_type pos = s.find_first_of('[');
+    string name;
+    PropertySet opts;
+    if( pos == string::npos ) {
+        name = s;
+    } else {
+        name = s.substr(0,pos);
 
-    // attach second part, merging properties if necessary
-    if( pos != string::npos ) {
-        if( s.at(s.length()-1) == ']' ) {
-            s = s.erase(s.length()-1,1) + ',' + _s.substr(pos+1);
-        } else
-            s = s + _s.substr(pos);
+        stringstream ss;
+        ss << s.substr(pos,s.length());
+        ss >> opts;
     }
+    return make_pair(name,opts);
+}
 
-    pair<string, PropertySet> result;
-    string & name = result.first;
-    PropertySet & opts = result.second;
 
-    pos = s.find_first_of('[');
-    if( pos == string::npos )
-        throw "Malformed method";
-    name = s.substr( 0, pos );
+pair<string, PropertySet> parseMethod( const string &_s, const map<string,string> & aliases ) {
+    // break string into method[properties]
+    pair<string,PropertySet> ps = parseMethodRaw(_s);
+    bool looped = false;
+
+    // as long as 'method' is an alias, update:
+    while( aliases.find(ps.first) != aliases.end() && !looped ) {
+        string astr = aliases.find(ps.first)->second;
+        pair<string,PropertySet> aps = parseMethodRaw(astr);
+        if( aps.first == ps.first )
+            looped = true;
+        // override aps properties by ps properties
+        aps.second.Set( ps.second );
+        // replace ps by aps
+        ps = aps;
+        // repeat until method name == alias name ('looped'), or
+        // there is no longer an alias 'method'
+    }
+    
+    // check whether name is valid
     size_t n = 0;
     for( ; strlen( DAINames[n] ) != 0; n++ )
-        if( name == DAINames[n] )
+        if( ps.first == DAINames[n] )
             break;
-    if( strlen( DAINames[n] ) == 0 && (name != "LDPC") )
-        DAI_THROW(UNKNOWN_DAI_ALGORITHM);
+    if( strlen( DAINames[n] ) == 0 && (ps.first != "LDPC") )
+        throw string("Unknown DAI algorithm \"") + ps.first + string("\" in \"") + _s + string("\"");
 
-    stringstream ss;
-    ss << s.substr(pos,s.length());
-    ss >> opts;
-    
-    return result;
+    return ps;
 }