Fixed testem failure caused by rounding error
[libdai.git] / include / dai / prob.h
index 4d6ea84..4b85374 100644 (file)
@@ -37,6 +37,7 @@
 #include <numeric>
 #include <functional>
 #include <dai/util.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
@@ -54,6 +55,11 @@ template <typename T> class TProb {
         std::vector<T> _p;
 
     public:
+        /// Iterator over entries
+       typedef typename std::vector<T>::iterator iterator;
+        /// Const iterator over entries
+       typedef typename std::vector<T>::const_iterator const_iterator;
+
         /// Enumerates different ways of normalizing a probability measure.
         /** 
          *  - NORMPROB means that the sum of all entries should be 1;
@@ -82,7 +88,7 @@ template <typename T> class TProb {
         /** \tparam Iterator Iterates over instances that can be cast to T.
          *  \param begin Points to first instance to be added.
          *  \param end Points just beyond last instance to be added.
-         *  \param sizeHint For efficiency, the number of elements can be speficied by sizeHint.
+         *  \param sizeHint For efficiency, the number of entries can be speficied by sizeHint.
          */
         template <typename Iterator>
         TProb( Iterator begin, Iterator end, size_t sizeHint=0 ) : _p() {
@@ -107,6 +113,18 @@ template <typename T> class TProb {
         
         /// Returns reference to the i'th entry
         T& operator[]( size_t i ) { return _p[i]; }
+        
+        /// Returns iterator pointing to first entry
+        iterator begin() { return _p.begin(); }
+
+        /// Returns const iterator pointing to first entry
+        const_iterator begin() const { return _p.begin(); }
+
+        /// Returns iterator pointing beyond last entry
+        iterator end() { return _p.end(); }
+
+        /// Returns const iterator pointing beyond last entry
+        const_iterator end() const { return _p.end(); }
 
         /// Sets all entries to x
         TProb<T> & fill(T x) { 
@@ -132,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 ) {
@@ -346,7 +370,7 @@ template <typename T> class TProb {
             TProb<T> x;
             x._p.reserve( size() );
             for( size_t i = 0; i < size(); i++ )
-                x._p.push_back( _p[i] < 0 ? (-p[i]) : p[i] );
+                x._p.push_back( _p[i] < 0 ? (-_p[i]) : _p[i] );
             return x;
         }
 
@@ -385,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;
@@ -402,28 +434,41 @@ 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 std::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();
-#ifdef DAI_DEBUG
-            assert( Z != 0.0 );
-#endif
-            *this /= Z;
+            if( Z == 0.0 )
+                DAI_THROW(NOT_NORMALIZABLE);
+            else
+                *this /= Z;
             return Z;
         }
 
@@ -436,7 +481,13 @@ template <typename T> class TProb {
     
         /// Returns true if one or more entries are NaN
         bool hasNaNs() const {
-            return (std::find_if( _p.begin(), _p.end(), isnan ) != _p.end());
+            bool foundnan = false;
+            for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
+                if( isnan( *x ) ) {
+                    foundnan = true;
+                    break;
+                }
+            return foundnan;
         }
 
         /// Returns true if one or more entries are negative
@@ -454,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];