Merged factor.h from SVN head

diff --git a/ChangeLog b/ChangeLog
index 00460a0..310044a 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -20,7 +20,9 @@ libDAI-0.2.2 (2008-??-??)
   - Replaced multind by Permute
   - Added MultiFor
   - Added State
+* Improved factor.h (merged from SVN head): mainly new functionality
 * Moved everything into namespace "dai"
+* Added ExactInf class for brute force exact inference
 * Renamed DEBUG to DAI_DEBUG to avoid conflicts
 * Added conditional compilation of inference methods
 * Contributions by Peter Gober:
@@ -93,7 +95,9 @@ libDAI-0.2.2 (2008-??-??)
       max() -> maxDiff()
       max_size() -> maxSize()
   - Prob::max() -> Prob::maxVal()
-  - Factor::max() -> Factor::maxVal()
+  - Factor::
+      max() -> maxVal()
+      part_sum() -> partSum()
   - toc() in util.h now returns seconds as a double
   - VarSet::operator&&
   - Properties -> PropertySet

diff --git a/STATUS b/STATUS
index 5b19fe5..97c59a1 100644 (file)
--- a/STATUS
+++ b/STATUS
@@ -17,8 +17,6 @@ A disadvantage of this approach may be worse cachability.
   to replace the linear searches that are performed every time for findVar()?
   No, a better idea is to avoid calls to findVar() as much as possible.
 - Can the FactorGraph constructors be optimized further?
-- Solve the proliferation of type checks for all different ENUM's in
-properties.cpp
 - Remove x2x?
 
 TODO FOR RELEASE:
@@ -54,14 +52,13 @@ alldai.h
 alldai.cpp
 properties.h
 properties.cpp
+factor.h
 
 FILES IN SVN HEAD THAT ARE STILL RELEVANT:
 ChangeLog
 README
 TODO
 prob.h
-factor.h
-factor.cpp
 factorgraph.h
 factorgraph.cpp
 regiongraph.h

diff --git a/include/dai/factor.h b/include/dai/factor.h
index 6f2e5b0..62b0943 100644 (file)
--- a/include/dai/factor.h
+++ b/include/dai/factor.h
@@ -41,6 +41,9 @@ typedef TFactor<Real>           Factor;
 // predefine friends
 template<typename T> Real           dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt );
 template<typename T> Real           KL_dist( const TFactor<T> & p, const TFactor<T> & q );
+template<typename T> Real           MutualInfo( const TFactor<T> & p );
+template<typename T> TFactor<T>     max( const TFactor<T> & P, const TFactor<T> & Q );
+template<typename T> TFactor<T>     min( const TFactor<T> & P, const TFactor<T> & Q );
 template<typename T> std::ostream&  operator<< (std::ostream& os, const TFactor<T>& P);
 
         
@@ -51,9 +54,9 @@ template <typename T> class TFactor {
         TProb<T>    _p;
 
     public:
-        // Default constructor
-        TFactor () : _vs(), _p(1,1.0) {}
-        
+        // Construct Factor with empty VarSet but nonempty _p
+        TFactor ( Real p = 1.0 ) : _vs(), _p(1,p) {}
+
         // Construct Factor from VarSet
         TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.states()) {}
         
@@ -64,7 +67,7 @@ template <typename T> class TFactor {
         TFactor( const VarSet& ns, const Real* p ) : _vs(ns), _p(_vs.states(),p) {}
 
         // Construct Factor from VarSet and TProb<T>
-        TFactor( const VarSet& ns, const TProb<T> p ) : _vs(ns), _p(p) {
+        TFactor( const VarSet& ns, const TProb<T>& p ) : _vs(ns), _p(p) {
 #ifdef DAI_DEBUG
             assert( _vs.states() == _p.size() );
 #endif
@@ -120,7 +123,9 @@ template <typename T> class TFactor {
             return *this;
         }
         TFactor<T> operator* (const TFactor<T>& Q) const;
+        TFactor<T> operator/ (const TFactor<T>& Q) const;
         TFactor<T>& operator*= (const TFactor<T>& Q) { return( *this = (*this * Q) ); }
+        TFactor<T>& operator/= (const TFactor<T>& Q) { return( *this = (*this / Q) ); }
         TFactor<T> operator+ (const TFactor<T>& Q) const {
 #ifdef DAI_DEBUG
             assert( Q._vs == _vs );
@@ -151,6 +156,24 @@ template <typename T> class TFactor {
             _p -= Q._p;
             return *this;
         }
+        TFactor<T>& operator+= (T q) { 
+            _p += q;
+            return *this;
+        }
+        TFactor<T>& operator-= (T q) { 
+            _p -= q;
+            return *this;
+        }
+        TFactor<T> operator+ (T q) const {
+            TFactor<T> result(*this); 
+            result._p += q; 
+            return result; 
+        }
+        TFactor<T> operator- (T q) const {
+            TFactor<T> result(*this); 
+            result._p -= q; 
+            return result; 
+        }
 
         TFactor<T> operator^ (Real a) const { TFactor<T> x; x._vs = _vs; x._p = _p^a; return x; }
         TFactor<T>& operator^= (Real a) { _p ^= a; return *this; }
@@ -159,6 +182,11 @@ template <typename T> class TFactor {
             _p.makeZero( epsilon );
             return *this;
         }
+
+        TFactor<T>& makePositive( Real epsilon ) {
+            _p.makePositive( epsilon );
+            return *this;
+        }
             
         TFactor<T> inverse() const { 
             TFactor<T> inv; 
@@ -191,6 +219,13 @@ template <typename T> class TFactor {
             return e; 
         }
 
+        TFactor<T> abs() const { 
+            TFactor<T> e; 
+            e._vs = _vs; 
+            e._p = _p.abs(); 
+            return e; 
+        }
+
         TFactor<T> log() const {
             TFactor<T> l; 
             l._vs = _vs; 
@@ -205,8 +240,8 @@ template <typename T> class TFactor {
             return l0; 
         }
 
-        T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); }
-        TFactor<T> normalized( typename Prob::NormType norm ) const { 
+        T normalize( typename Prob::NormType norm = Prob::NORMPROB ) { return _p.normalize( norm ); }
+        TFactor<T> normalized( typename Prob::NormType norm = Prob::NORMPROB ) const { 
             TFactor<T> result;
             result._vs = _vs;
             result._p = _p.normalized( norm );
@@ -229,16 +264,29 @@ template <typename T> class TFactor {
             return result;
         }
 
-        // returns unnormalized marginal
-        TFactor<T> part_sum(const VarSet & ns) const;
-        // returns normalized marginal
-        TFactor<T> marginal(const VarSet & ns) const { return part_sum(ns).normalized( Prob::NORMPROB ); }
+        // returns unnormalized marginal; ns should be a subset of vars()
+        TFactor<T> partSum(const VarSet & ns) const;
+        // returns (normalized by default) marginal; ns should be a subset of vars()
+        TFactor<T> marginal(const VarSet & ns, bool normed = true) const { if(normed) return partSum(ns).normalized(); else return partSum(ns); }
+        // sums out all variables except those in ns
+        TFactor<T> notSum(const VarSet & ns) const { return partSum(vars() ^ ns); }
+
+        // embeds this factor in larger varset ns
+        TFactor<T> embed(const VarSet & ns) const { 
+            VarSet vs = vars();
+            assert( ns >> vs );
+            if( vs == ns )
+                return *this;
+            else
+                return (*this) * Factor(ns / vs, 1.0);
+        }
 
         bool hasNaNs() const { return _p.hasNaNs(); }
         bool hasNegatives() const { return _p.hasNegatives(); }
         T totalSum() const { return _p.totalSum(); }
         T maxAbs() const { return _p.maxAbs(); }
         T maxVal() const { return _p.maxVal(); }
+        T minVal() const { return _p.minVal(); }
         Real entropy() const { return _p.entropy(); }
         T strength( const Var &i, const Var &j ) const;
 
@@ -253,11 +301,12 @@ template <typename T> class TFactor {
             }
         }
         friend Real KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
+        friend Real MutualInfo <> ( const TFactor<T> & P );
         template<class U> friend std::ostream& operator<< (std::ostream& os, const TFactor<U>& P);
 };
 
 
-template<typename T> TFactor<T> TFactor<T>::part_sum(const VarSet & ns) const {
+template<typename T> TFactor<T> TFactor<T>::partSum(const VarSet & ns) const {
 #ifdef DAI_DEBUG
     assert( ns << _vs );
 #endif
@@ -294,6 +343,19 @@ template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& Q) cons
 }
 
 
+template<typename T> TFactor<T> TFactor<T>::operator/ (const TFactor<T>& Q) const {
+    TFactor<T> quot( _vs + Q._vs, 0.0 );
+
+    IndexFor i1(_vs, quot._vs);
+    IndexFor i2(Q._vs, quot._vs);
+
+    for( size_t i = 0; i < quot._p.size(); i++, ++i1, ++i2 )
+        quot._p[i] += _p[i1] / Q._p[i2];
+
+    return quot;
+}
+
+
 template<typename T> Real KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
     if( P._vs.empty() || Q._vs.empty() )
         return -1;
@@ -306,6 +368,26 @@ template<typename T> Real KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
 }
 
 
+// calculate mutual information of x_i and x_j where P.vars() = \{x_i,x_j\}
+template<typename T> Real MutualInfo(const TFactor<T> & P) {
+    assert( P._vs.size() == 2 );
+    VarSet::const_iterator it = P._vs.begin();
+    Var i = *it; it++; Var j = *it;
+    TFactor<T> projection = P.marginal(i) * P.marginal(j);
+    return real( KL_dist( P.normalized(), projection ) );
+}
+
+
+template<typename T> TFactor<T> max( const TFactor<T> & P, const TFactor<T> & Q ) {
+    assert( P._vs == Q._vs );
+    return TFactor<T>( P._vs, min( P.p(), Q.p() ) );
+}
+
+template<typename T> TFactor<T> min( const TFactor<T> & P, const TFactor<T> & Q ) {
+    assert( P._vs == Q._vs );
+    return TFactor<T>( P._vs, max( P.p(), Q.p() ) );
+}
+
 // calculate N(psi, i, j)
 template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
 #ifdef DAI_DEBUG
@@ -313,7 +395,7 @@ template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const
     assert( _vs.contains( j ) );
     assert( i != j );
 #endif
-    VarSet ij = i | j;
+    VarSet ij(i, j);
 
     T max = 0.0;
     for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
@@ -344,8 +426,8 @@ template<typename T> TFactor<T> RemoveFirstOrderInteractions( const TFactor<T> &
     VarSet vars = psi.vars();
     for( size_t iter = 0; iter < 100; iter++ ) {
         for( VarSet::const_iterator n = vars.begin(); n != vars.end(); n++ )
-            result = result * result.part_sum(*n).inverse();
-        result.normalize( Prob::NORMPROB );
+            result = result * result.partSum(*n).inverse();
+        result.normalize();
     }
 
     return result;

diff --git a/src/jtree.cpp b/src/jtree.cpp
index 558e367..7de6a36 100644 (file)
--- a/src/jtree.cpp
+++ b/src/jtree.cpp
@@ -215,7 +215,7 @@ void JTree::runHUGIN() {
     for( size_t i = _RTree.size(); (i--) != 0; ) {
 //      Make outer region _RTree[i].n1 consistent with outer region _RTree[i].n2
 //      IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
-        Factor new_Qb = _Qa[_RTree[i].n2].part_sum( IR( i ) );
+        Factor new_Qb = _Qa[_RTree[i].n2].partSum( IR( i ) );
         _logZ += log(new_Qb.normalize( Prob::NORMPROB ));
         _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] ); 
         _Qb[i] = new_Qb;
@@ -256,7 +256,7 @@ void JTree::runShaferShenoy() {
         foreach( const Neighbor &k, nbOR(i) )
             if( k != e ) 
                 piet *= message( i, k.iter );
-        message( j, _e ) = piet.part_sum( IR(e) );
+        message( j, _e ) = piet.partSum( IR(e) );
         _logZ += log( message(j,_e).normalize( Prob::NORMPROB ) );
     }
 
@@ -494,7 +494,7 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
                             _Qa[T[i].n2] *= piet; 
                         }
 
-                    Factor new_Qb = _Qa[T[i].n2].part_sum( IR( b[i] ) );
+                    Factor new_Qb = _Qa[T[i].n2].partSum( IR( b[i] ) );
                     logZ += log(new_Qb.normalize( Prob::NORMPROB ));
                     _Qa[T[i].n1] *= new_Qb.divided_by( _Qb[b[i]] ); 
                     _Qb[b[i]] = new_Qb;
@@ -503,7 +503,7 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
 
                 Factor piet( nsrem, 0.0 );
                 piet[s] = exp(logZ);
-                Pns += piet * _Qa[T[0].n1].part_sum( ns / nsrem );      // OPTIMIZE ME
+                Pns += piet * _Qa[T[0].n1].partSum( ns / nsrem );      // OPTIMIZE ME
 
                 // Restore clamped beliefs
                 for( map<size_t,Factor>::const_iterator alpha = _Qa_old.begin(); alpha != _Qa_old.end(); alpha++ )

diff --git a/src/lc.cpp b/src/lc.cpp
index f13dead..7e6996c 100644 (file)
--- a/src/lc.cpp
+++ b/src/lc.cpp
@@ -258,10 +258,10 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
     Factor A_I;
     for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
         if( var(i) != *k )
-            A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).part_sum( Ivars / var(i) );
+            A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).partSum( Ivars / var(i) );
     if( Ivars.size() > 1 )
         A_I ^= (1.0 / (Ivars.size() - 1));
-    Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).part_sum( Ivars / var(i) );
+    Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).partSum( Ivars / var(i) );
     Factor quot = A_I.divided_by(A_Ii);
 
     piet *= quot.divided_by( _phis[i][_I] ).normalized( Prob::NORMPROB );

diff --git a/src/mf.cpp b/src/mf.cpp
index 582bd57..526c776 100644 (file)
--- a/src/mf.cpp
+++ b/src/mf.cpp
@@ -112,7 +112,7 @@ double MF::run() {
                     henk *= _beliefs[j];
             piet = factor(I).log0();
             piet *= henk;
-            piet = piet.part_sum(var(i));
+            piet = piet.partSum(var(i));
             piet = piet.exp();
             jan *= piet; 
         }

diff --git a/src/treeep.cpp b/src/treeep.cpp
index a44007f..8be99fe 100644 (file)
--- a/src/treeep.cpp
+++ b/src/treeep.cpp
@@ -151,14 +151,14 @@ void TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &
                     delta[s(*n)] = 1.0;
                     _Qa[_RTree[i].n2] *= delta;
                 }
-            Factor new_Qb = _Qa[_RTree[i].n2].part_sum( _Qb[i].vars() );
+            Factor new_Qb = _Qa[_RTree[i].n2].partSum( _Qb[i].vars() );
             _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] ); 
             _Qb[i] = new_Qb;
         }
 
         // DistributeEvidence
         for( size_t i = 0; i < _RTree.size(); i++ ) {
-            Factor new_Qb = _Qa[_RTree[i].n1].part_sum( _Qb[i].vars() );
+            Factor new_Qb = _Qa[_RTree[i].n1].partSum( _Qb[i].vars() );
             _Qa[_RTree[i].n2] *= new_Qb.divided_by( _Qb[i] ); 
             _Qb[i] = new_Qb;
         }