Various cleanups and documentation improvements for SmallSet and Prob
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 19 Oct 2009 10:57:02 +0000 (12:57 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 19 Oct 2009 10:57:02 +0000 (12:57 +0200)
include/dai/factor.h
include/dai/hak.h
include/dai/jtree.h
include/dai/prob.h
include/dai/smallset.h
include/dai/util.h
src/hak.cpp
src/util.cpp
utils/createfg.cpp

index e97ea78..cdb324c 100644 (file)
@@ -283,19 +283,6 @@ template <typename T> class TFactor {
             return pointwiseOp(*this,f,std::minus<T>());
         }
 
-
-        /// Sets all values that are smaller than epsilon to 0
-        TFactor<T>& makeZero( T epsilon ) {
-            _p.makeZero( epsilon );
-            return *this;
-        }
-
-        /// Sets all values that are smaller than epsilon to epsilon
-        TFactor<T>& makePositive( T epsilon ) {
-            _p.makePositive( epsilon );
-            return *this;
-        }
-
         /// Returns pointwise inverse of *this.
         /** If zero == true, uses 1 / 0 == 0; otherwise 1 / 0 == Inf.
          */
index ccca19b..96a2f00 100644 (file)
@@ -46,10 +46,10 @@ class HAK : public DAIAlgRG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of possible cluster choices
-            DAI_ENUM(ClustersType,MIN,DELTA,LOOP)
+            DAI_ENUM(ClustersType,MIN,DELTA,LOOP);
 
             /// Enumeration of possible message initializations
-            DAI_ENUM(InitType,UNIFORM,RANDOM)
+            DAI_ENUM(InitType,UNIFORM,RANDOM);
 
             /// Verbosity
             size_t verbose;
index 6762d34..8107757 100644 (file)
@@ -52,7 +52,7 @@ class JTree : public DAIAlgRG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of possible JTree updates
-            DAI_ENUM(UpdateType,HUGIN,SHSH)
+            DAI_ENUM(UpdateType,HUGIN,SHSH);
 
             /// Enumeration of inference variants
             DAI_ENUM(InfType,SUMPROD,MAXPROD);
index 302a4d6..b5ad484 100644 (file)
@@ -10,9 +10,8 @@
 
 
 /// \file
-/// \brief Defines TProb<T> and Prob classes
+/// \brief Defines TProb<T> and Prob classes which represent (probability) vectors
 /// \todo Rename to Vector<T>
-/// \todo Check documentation
 
 
 #ifndef __defined_libdai_prob_h
@@ -33,10 +32,13 @@ namespace dai {
 
 
 /// Represents a vector with entries of type \a T.
-/** A TProb<T> is a std::vector<T> with an interface designed for dealing with probability mass functions.
- *  It is mainly used for representing measures on a finite outcome space, e.g., the probability
- *  distribution of a discrete random variable.
- *  \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
+/** It is simply a <tt>std::vector</tt><<em>T</em>> with an interface designed for dealing with probability mass functions.
+ *
+ *  It is mainly used for representing measures on a finite outcome space, for example, the probability
+ *  distribution of a discrete random variable. However, entries are not necessarily non-negative; it is also used to
+ *  represent logarithms of probability mass functions.
+ *
+ *  \tparam T Should be a scalar that is castable from and to dai::Real and should support elementary arithmetic operations.
  */
 template <typename T> class TProb {
     private:
@@ -44,11 +46,6 @@ 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;
@@ -57,36 +54,38 @@ template <typename T> class TProb {
         typedef enum { NORMPROB, NORMLINF } NormType;
         /// Enumerates different distance measures between probability measures.
         /**
-         *  - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
-         *  - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
-         *  - DISTTV is the Total Variation distance;
-         *  - DISTKL is the Kullback-Leibler distance.
+         *  - DISTL1 is the \f$\ell_1\f$ distance (sum of absolute values of pointwise difference);
+         *  - DISTLINF is the \f$\ell_\infty\f$ distance (maximum absolute value of pointwise difference);
+         *  - DISTTV is the total variation distance (half of the \f$\ell_1\f$ distance);
+         *  - DISTKL is the Kullback-Leibler distance (\f$\sum_i p_i (\log p_i - \log q_i)\f$).
          */
         typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
 
-        /// Default constructor
+    /// \name Constructors and destructors
+    //@{
+        /// Default constructor (constructs empty vector)
         TProb() : _p() {}
 
-        /// Construct uniform distribution over n outcomes, i.e., a vector of length n with each entry set to 1/n
-        explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
+        /// Construct uniform probability distribution over \a n outcomes (i.e., a vector of length \a n with each entry set to \f$1/n\f$)
+        explicit TProb( size_t n ) : _p(std::vector<T>(n, (T)1 / n)) {}
 
-        /// Construct vector of length n with each entry set to p
-        explicit TProb( size_t n, Real p ) : _p(n, (T)p) {}
+        /// Construct vector of length \a n with each entry set to \a p
+        explicit TProb( size_t n, T p ) : _p(n, p) {}
 
         /// Construct vector from a range
-        /** \tparam Iterator Iterates over instances that can be cast to T.
+        /** \tparam TIterator Iterates over instances that can be cast to \a 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 entries can be speficied by sizeHint.
+         *  \param sizeHint For efficiency, the number of entries can be speficied by \a sizeHint.
          */
-        template <typename Iterator>
-        TProb( Iterator begin, Iterator end, size_t sizeHint=0 ) : _p() {
+        template <typename TIterator>
+        TProb( TIterator begin, TIterator end, size_t sizeHint=0 ) : _p() {
             _p.reserve( sizeHint );
             _p.insert( _p.begin(), begin, end );
         }
 
-        /// Construct vector from a vector
-        /** \tparam S type of elements in v
+        /// Construct vector from another vector
+        /** \tparam S type of elements in \a v (should be castable to type \a T)
          *  \param v vector used for initialization
          */
         template <typename S>
@@ -94,14 +93,49 @@ template <typename T> class TProb {
             _p.reserve( v.size() );
             _p.insert( _p.begin(), v.begin(), v.end() );
         }
+    //@}
+
+        /// Constant iterator over the elements
+        typedef typename std::vector<T>::const_iterator const_iterator;
+        /// Iterator over the elements
+        typedef typename std::vector<T>::iterator iterator;
+        /// Constant reverse iterator over the elements
+        typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
+        /// Reverse iterator over the elements
+        typedef typename std::vector<T>::reverse_iterator reverse_iterator;
+
+    /// @name Iterator interface
+    //@{
+        /// Returns iterator that points to the first element
+        iterator begin() { return _p.begin(); }
+        /// Returns constant iterator that points to the first element
+        const_iterator begin() const { return _p.begin(); }
+
+        /// Returns iterator that points beyond the last element
+        iterator end() { return _p.end(); }
+        /// Returns constant iterator that points beyond the last element
+        const_iterator end() const { return _p.end(); }
+
+        /// Returns reverse iterator that points to the last element
+        reverse_iterator rbegin() { return _p.rbegin(); }
+        /// Returns constant reverse iterator that points to the last element
+        const_reverse_iterator rbegin() const { return _p.rbegin(); }
 
-        /// Returns a const reference to the vector
+        /// Returns reverse iterator that points beyond the first element
+        reverse_iterator rend() { return _p.rend(); }
+        /// Returns constant reverse iterator that points beyond the first element
+        const_reverse_iterator rend() const { return _p.rend(); }
+    //@}
+
+    /// \name Queries
+    //@{
+        /// Returns a const reference to the wrapped vector
         const std::vector<T> & p() const { return _p; }
 
-        /// Returns a reference to the vector
+        /// Returns a reference to the wrapped vector
         std::vector<T> & p() { return _p; }
 
-        /// Returns a copy of the i'th entry
+        /// Returns a copy of the \a i 'th entry
         T operator[]( size_t i ) const {
 #ifdef DAI_DEBUG
             return _p.at(i);
@@ -110,114 +144,101 @@ template <typename T> class TProb {
 #endif
         }
 
-        /// Returns reference to the i'th entry
+        /// Returns reference to the \a 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(); }
+        /// Returns length of the vector (i.e., the number of entries)
+        size_t size() const { return _p.size(); }
 
-        /// Sets all entries to x
-        TProb<T> & fill(T x) {
-            std::fill( _p.begin(), _p.end(), x );
-            return *this;
-        }
-
-        /// Draws all entries i.i.d. from a uniform distribution on [0,1)
-        TProb<T> & randomize() {
-            std::generate(_p.begin(), _p.end(), rnd_uniform);
-            return *this;
-        }
-
-        /// Returns length of the vector, i.e., the number of entries
-        size_t size() const {
-            return _p.size();
-        }
-
-        /// Sets entries that are smaller than epsilon to 0
-        TProb<T>& makeZero( Real epsilon ) {
+        /// Returns the Shannon entropy of *this, \f$-\sum_i p_i \log p_i\f$
+        T entropy() const {
+            T S = 0;
             for( size_t i = 0; i < size(); i++ )
-                if( fabs(_p[i]) < epsilon )
-                    _p[i] = 0;
-            return *this;
-        }
-
-        /// Set all entries to 1.0/size()
-        TProb<T>& setUniform () {
-            fill(1.0/size());
-            return *this;
+                S -= (_p[i] == 0 ? 0 : _p[i] * dai::log(_p[i]));
+            return S;
         }
 
-        /// Sets entries that are smaller than epsilon to epsilon
-        TProb<T>& makePositive( Real epsilon ) {
-            for( size_t i = 0; i < size(); i++ )
-                if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
-                    _p[i] = epsilon;
-            return *this;
+        /// Returns maximum value of all entries
+        T max() const {
+            T Z = *std::max_element( _p.begin(), _p.end() );
+            return Z;
         }
 
-        /// Multiplies each entry with scalar x
-        TProb<T>& operator*= (T x) {
-            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
-            return *this;
+        /// Returns minimum value of all entries
+        T min() const {
+            T Z = *std::min_element( _p.begin(), _p.end() );
+            return Z;
         }
 
-        /// Returns product of *this with scalar x
-        TProb<T> operator* (T x) const {
-            TProb<T> prod( *this );
-            prod *= x;
-            return prod;
+        /// Returns a pair consisting of the index of the maximum value and the maximum value itself
+        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);
         }
 
-        /// Divides each entry by scalar x
-        TProb<T>& operator/= (T x) {
-            DAI_DEBASSERT( x != 0.0 );
-            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
-            return *this;
+        /// Returns sum of all entries
+        T sum() const {
+            T Z = std::accumulate( _p.begin(),  _p.end(), (T)0 );
+            return Z;
         }
 
-        /// Returns quotient of *this and scalar x
-        TProb<T> operator/ (T x) const {
-            TProb<T> quot( *this );
-            quot /= x;
-            return quot;
+        /// Return sum of absolute value of all entries
+        T sumAbs() const {
+            T s = 0;
+            for( size_t i = 0; i < size(); i++ )
+                s += dai::abs(_p[i]);
+            return s;
         }
 
-        /// Adds scalar x to each entry
-        TProb<T>& operator+= (T x) {
-            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
-            return *this;
+        /// Returns maximum absolute value of all entries
+        T maxAbs() const {
+            T Z = 0;
+            for( size_t i = 0; i < size(); i++ ) {
+                T mag = dai::abs(_p[i]);
+                if( mag > Z )
+                    Z = mag;
+            }
+            return Z;
         }
 
-        /// Returns sum of *this and scalar x
-        TProb<T> operator+ (T x) const {
-            TProb<T> sum( *this );
-            sum += x;
-            return sum;
+        /// Returns true if one or more entries are NaN
+        bool hasNaNs() const {
+            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;
         }
 
-        /// Subtracts scalar x from each entry
-        TProb<T>& operator-= (T x) {
-            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
-            return *this;
+        /// Returns true if one or more entries are negative
+        bool hasNegatives() const {
+            return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<T>(), (T)0 ) ) != _p.end());
         }
 
-        /// Returns difference of *this and scalar x
-        TProb<T> operator- (T x) const {
-            TProb<T> diff( *this );
-            diff -= x;
-            return diff;
+        /// Returns a random index, according to the (normalized) distribution described by *this
+        size_t draw() {
+            Real x = rnd_uniform() * sum();
+            T s = 0;
+            for( size_t i = 0; i < size(); i++ ) {
+                s += _p[i];
+                if( s > x )
+                    return i;
+            }
+            return( size() - 1 );
         }
 
-        /// Lexicographical comparison (sizes should be identical)
+        /// Lexicographical comparison
+        /** \pre this->size() == q.size()
+         */
         bool operator<= (const TProb<T> & q) const {
             DAI_DEBASSERT( size() == q.size() );
             for( size_t i = 0; i < size(); i++ )
@@ -225,108 +246,10 @@ template <typename T> class TProb {
                     return false;
             return true;
         }
+    //@}
 
-        /// Pointwise multiplication with q (sizes should be identical)
-        TProb<T>& operator*= (const TProb<T> & q) {
-            DAI_DEBASSERT( size() == q.size() );
-            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
-            return *this;
-        }
-
-        /// Return product of *this with q (sizes should be identical)
-        TProb<T> operator* (const TProb<T> & q) const {
-            DAI_DEBASSERT( size() == q.size() );
-            TProb<T> prod( *this );
-            prod *= q;
-            return prod;
-        }
-
-        /// Pointwise addition with q (sizes should be identical)
-        TProb<T>& operator+= (const TProb<T> & q) {
-            DAI_DEBASSERT( size() == q.size() );
-            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
-            return *this;
-        }
-
-        /// Returns sum of *this and q (sizes should be identical)
-        TProb<T> operator+ (const TProb<T> & q) const {
-            DAI_DEBASSERT( size() == q.size() );
-            TProb<T> sum( *this );
-            sum += q;
-            return sum;
-        }
-
-        /// Pointwise subtraction of q (sizes should be identical)
-        TProb<T>& operator-= (const TProb<T> & q) {
-            DAI_DEBASSERT( size() == q.size() );
-            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
-            return *this;
-        }
-
-        /// Return *this minus q (sizes should be identical)
-        TProb<T> operator- (const TProb<T> & q) const {
-            DAI_DEBASSERT( size() == q.size() );
-            TProb<T> diff( *this );
-            diff -= q;
-            return diff;
-        }
-
-        /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
-        TProb<T>& operator/= (const TProb<T> & q) {
-            DAI_DEBASSERT( size() == q.size() );
-            for( size_t i = 0; i < size(); i++ ) {
-                if( q[i] == 0.0 )
-                    _p[i] = 0.0;
-                else
-                    _p[i] /= q[i];
-            }
-            return *this;
-        }
-
-        /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
-        TProb<T>& divide (const TProb<T> & q) {
-            DAI_DEBASSERT( size() == q.size() );
-            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
-            return *this;
-        }
-
-        /// Returns quotient of *this with q (sizes should be identical)
-        TProb<T> operator/ (const TProb<T> & q) const {
-            DAI_DEBASSERT( size() == q.size() );
-            TProb<T> quot( *this );
-            quot /= q;
-            return quot;
-        }
-
-        /// Returns pointwise inverse
-        /** If zero==true; uses 1/0==0, otherwise 1/0==Inf.
-         */
-        TProb<T> inverse(bool zero=true) const {
-            TProb<T> inv;
-            inv._p.reserve( size() );
-            if( zero )
-                for( size_t i = 0; i < size(); i++ )
-                    inv._p.push_back( _p[i] == 0.0 ? 0.0 : 1.0 / _p[i] );
-            else
-                for( size_t i = 0; i < size(); i++ )
-                    inv._p.push_back( 1.0 / _p[i] );
-            return inv;
-        }
-
-        /// Raises entries to the power a
-        TProb<T>& operator^= (Real a) {
-            if( a != 1.0 )
-                std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::ptr_fun<T, Real, T>(std::pow), a) );
-            return *this;
-        }
-
-        /// Returns *this raised to the power a
-        TProb<T> operator^ (Real a) const {
-            TProb<T> power(*this);
-            power ^= a;
-            return power;
-        }
-
+    /// \name Unary transformations
+    //@{
         /// Returns pointwise signum
         TProb<T> sgn() const {
             TProb<T> x;
@@ -347,36 +270,18 @@ 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( abs(_p[i]) );
             return x;
         }
 
-        /// Applies exp pointwise
-        const TProb<T>& takeExp() {
-            std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::exp) );
-            return *this;
-        }
-
-        /// Applies log pointwise
-        /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
-         */
-        const TProb<T>& takeLog(bool zero=false) {
-            if( zero ) {
-                for( size_t i = 0; i < size(); i++ )
-                   _p[i] = ( (_p[i] == 0.0) ? 0.0 : std::log( _p[i] ) );
-            } else
-                std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::log) );
-            return *this;
-        }
-
-        /// Returns pointwise exp
+        /// Returns pointwise exponent
         TProb<T> exp() const {
             TProb<T> e(*this);
             e.takeExp();
             return e;
         }
 
-        /// Returns pointwise log
+        /// Returns pointwise logarithm
         /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
          */
         TProb<T> log(bool zero=false) const {
@@ -385,54 +290,61 @@ template <typename T> class TProb {
             return l;
         }
 
-        /// Returns sum of all entries
-        T sum() const {
-            T Z = std::accumulate( _p.begin(),  _p.end(), (T)0 );
-            return Z;
+        /// Returns pointwise inverse
+        /** If zero==true; uses 1/0==0, otherwise 1/0==Inf.
+         */
+        TProb<T> inverse(bool zero=true) const {
+            TProb<T> inv;
+            inv._p.reserve( size() );
+            if( zero )
+                for( size_t i = 0; i < size(); i++ )
+                    inv._p.push_back( _p[i] == 0.0 ? 0.0 : 1.0 / _p[i] );
+            else
+                for( size_t i = 0; i < size(); i++ )
+                    inv._p.push_back( 1.0 / _p[i] );
+            return inv;
         }
 
-        /// 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 normalized copy of \c *this, using the specified norm
+        TProb<T> normalized( NormType norm = NORMPROB ) const {
+            TProb<T> result(*this);
+            result.normalize( norm );
+            return result;
         }
+    //@}
 
-        /// Returns maximum absolute value of all entries
-        T maxAbs() const {
-            T Z = 0;
-            for( size_t i = 0; i < size(); i++ ) {
-                Real mag = fabs( (Real) _p[i] );
-                if( mag > Z )
-                    Z = mag;
-            }
-            return Z;
+    /// \name Unary operations
+    //@{
+        /// Draws all entries i.i.d. from a uniform distribution on [0,1)
+        TProb<T>& randomize() {
+            std::generate( _p.begin(), _p.end(), rnd_uniform );
+            return *this;
         }
 
-        /// Returns maximum value of all entries
-        T max() const {
-            T Z = *std::max_element( _p.begin(), _p.end() );
-            return Z;
+        /// Sets all entries to \f$1/n\f$ where \a n is the length of the vector
+        TProb<T>& setUniform () {
+            fill( (T)1 / size() );
+            return *this;
         }
 
-        /// Returns minimum value of all entries
-        T min() const {
-            T Z = *std::min_element( _p.begin(), _p.end() );
-            return Z;
+        /// Applies exponent pointwise
+        const TProb<T>& takeExp() {
+            for( size_t i = 0; i < size(); i++ )
+               _p[i] = dai::exp(_p[i]);
+            return *this;
         }
 
-        /// 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);
+        /// Applies logarithm pointwise
+        /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
+         */
+        const TProb<T>& takeLog(bool zero=false) {
+            if( zero ) {
+                for( size_t i = 0; i < size(); i++ )
+                   _p[i] = ( (_p[i] == 0.0) ? 0.0 : dai::log(_p[i]) );
+            } else
+                for( size_t i = 0; i < size(); i++ )
+                   _p[i] = dai::log(_p[i]);
+            return *this;
         }
 
         /// Normalizes vector using the specified norm
@@ -448,67 +360,215 @@ template <typename T> class TProb {
                 *this /= Z;
             return Z;
         }
+    //@}
 
-        /// Returns normalized copy of *this, using the specified norm
-        TProb<T> normalized( NormType norm = NORMPROB ) const {
-            TProb<T> result(*this);
-            result.normalize( norm );
-            return result;
+    /// \name Operations with scalars
+    //@{
+        /// Sets all entries to \a x
+        TProb<T> & fill(T x) {
+            std::fill( _p.begin(), _p.end(), x );
+            return *this;
         }
 
-        /// Returns true if one or more entries are NaN
-        bool hasNaNs() const {
-            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;
+        /// Adds scalar \a x to each entry
+        TProb<T>& operator+= (T x) {
+            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
+            return *this;
         }
 
-        /// Returns true if one or more entries are negative
-        bool hasNegatives() const {
-            return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
+        /// Subtracts scalar \a x from each entry
+        TProb<T>& operator-= (T x) {
+            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
+            return *this;
         }
 
-        /// Returns entropy of *this
-        Real entropy() const {
-            Real S = 0.0;
-            for( size_t i = 0; i < size(); i++ )
-                S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));
-            return S;
+        /// Multiplies each entry with scalar \a x
+        TProb<T>& operator*= (T x) {
+            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
+            return *this;
         }
 
-        /// Returns a random index, according to the (normalized) distribution described by *this
-        size_t draw() {
-            double x = rnd_uniform() * sum();
-            T s = 0;
+        /// Divides each entry by scalar \a x
+        TProb<T>& operator/= (T x) {
+            DAI_DEBASSERT( x != 0 );
+            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
+            return *this;
+        }
+
+        /// Raises entries to the power \a x
+        TProb<T>& operator^= (T x) {
+            if( x != (T)1 )
+                std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::ptr_fun<T, T, T>(std::pow), x) );
+            return *this;
+        }
+    //@}
+
+    /// \name Transformations with scalars
+    //@{
+        /// Returns sum of \c *this and scalar \a x
+        TProb<T> operator+ (T x) const {
+            TProb<T> sum( *this );
+            sum += x;
+            return sum;
+        }
+
+        /// Returns difference of \c *this and scalar \a x
+        TProb<T> operator- (T x) const {
+            TProb<T> diff( *this );
+            diff -= x;
+            return diff;
+        }
+
+        /// Returns product of \c *this with scalar \a x
+        TProb<T> operator* (T x) const {
+            TProb<T> prod( *this );
+            prod *= x;
+            return prod;
+        }
+
+        /// Returns quotient of \c *this and scalar \a x
+        TProb<T> operator/ (T x) const {
+            TProb<T> quot( *this );
+            quot /= x;
+            return quot;
+        }
+
+        /// Returns *this raised to the power \a x
+        TProb<T> operator^ (T x) const {
+            TProb<T> power(*this);
+            power ^= x;
+            return power;
+        }
+    //@}
+
+    /// \name Operations with other equally-sized vectors
+    //@{
+        /// Pointwise addition with \a q
+        /** \pre this->size() == q.size()
+         */
+        TProb<T>& operator+= (const TProb<T> & q) {
+            DAI_DEBASSERT( size() == q.size() );
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
+            return *this;
+        }
+
+        /// Pointwise subtraction of \a q
+        /** \pre this->size() == q.size()
+         */
+        TProb<T>& operator-= (const TProb<T> & q) {
+            DAI_DEBASSERT( size() == q.size() );
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
+            return *this;
+        }
+
+        /// Pointwise multiplication with \a q
+        /** \pre this->size() == q.size()
+         */
+        TProb<T>& operator*= (const TProb<T> & q) {
+            DAI_DEBASSERT( size() == q.size() );
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
+            return *this;
+        }
+
+        /// Pointwise division by \a q, where division by 0 yields 0
+        /** \pre this->size() == q.size()
+         *  \see divide(const TProb<T> &)
+         */
+        TProb<T>& operator/= (const TProb<T> & q) {
+            DAI_DEBASSERT( size() == q.size() );
             for( size_t i = 0; i < size(); i++ ) {
-                s += _p[i];
-                if( s > x )
-                    return i;
+                if( q[i] == 0.0 )
+                    _p[i] = 0.0;
+                else
+                    _p[i] /= q[i];
             }
-            return( size() - 1 );
+            return *this;
+        }
+
+        /// Pointwise division by \a q, where division by 0 yields +Inf
+        /** \pre this->size() == q.size()
+         *  \see operator/=(const TProb<T> &)
+         */
+        TProb<T>& divide (const TProb<T> & q) {
+            DAI_DEBASSERT( size() == q.size() );
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
+            return *this;
+        }
+    //@}
+
+    /// \name Transformations with other equally-sized vectors
+    //@{
+        /// Returns sum of \c *this and \a q
+        /** \pre this->size() == q.size()
+         */
+        TProb<T> operator+ (const TProb<T> & q) const {
+            DAI_DEBASSERT( size() == q.size() );
+            TProb<T> sum( *this );
+            sum += q;
+            return sum;
+        }
+
+        /// Return \c *this minus \a q
+        /** \pre this->size() == q.size()
+         */
+        TProb<T> operator- (const TProb<T> & q) const {
+            DAI_DEBASSERT( size() == q.size() );
+            TProb<T> diff( *this );
+            diff -= q;
+            return diff;
+        }
+
+        /// Return product of \c *this with \a q
+        /** \pre this->size() == q.size()
+         */
+        TProb<T> operator* (const TProb<T> & q) const {
+            DAI_DEBASSERT( size() == q.size() );
+            TProb<T> prod( *this );
+            prod *= q;
+            return prod;
+        }
+
+        /// Returns quotient of \c *this with \a q, where division by 0 yields +Inf
+        /** \pre this->size() == q.size()
+         *  \see divided_by(const TProb<T> &)
+         */
+        TProb<T> operator/ (const TProb<T> & q) const {
+            DAI_DEBASSERT( size() == q.size() );
+            TProb<T> quot( *this );
+            quot /= q;
+            return quot;
+        }
+
+        /// Pointwise division by \a q, where division by 0 yields 0
+        /** \pre this->size() == q.size()
+         *  \see operator/(const TProb<T> &)
+         */
+        TProb<T> divided_by (const TProb<T> & q) const {
+            DAI_DEBASSERT( size() == q.size() );
+            TProb<T> quot( *this );
+            quot.divide(q);
+            return quot;
         }
+    //@}
 };
 
 
-/// Returns distance of p and q (sizes should be identical), measured using distance measure dt
+/// Returns distance between \a p and \a q, measured using distance measure \a dt
 /** \relates TProb
+ *  \pre this->size() == q.size()
  */
-template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
+template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
     DAI_DEBASSERT( p.size() == q.size() );
-    Real result = 0.0;
+    T result = 0;
     switch( dt ) {
         case TProb<T>::DISTL1:
             for( size_t i = 0; i < p.size(); i++ )
-                result += fabs((Real)p[i] - (Real)q[i]);
+                result += abs(p[i] - q[i]);
             break;
 
         case TProb<T>::DISTLINF:
             for( size_t i = 0; i < p.size(); i++ ) {
-                Real z = fabs((Real)p[i] - (Real)q[i]);
+                T z = abs(p[i] - q[i]);
                 if( z > result )
                     result = z;
             }
@@ -516,14 +576,14 @@ template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename T
 
         case TProb<T>::DISTTV:
             for( size_t i = 0; i < p.size(); i++ )
-                result += fabs((Real)p[i] - (Real)q[i]);
-            result *= 0.5;
+                result += abs(p[i] - q[i]);
+            result /= 2;
             break;
 
         case TProb<T>::DISTKL:
             for( size_t i = 0; i < p.size(); i++ ) {
                 if( p[i] != 0.0 )
-                    result += p[i] * (std::log(p[i]) - std::log(q[i]));
+                    result += p[i] * (dai::log(p[i]) - dai::log(q[i]));
             }
     }
     return result;
@@ -571,7 +631,7 @@ template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
 }
 
 
-/// Represents a vector with entries of type Real.
+/// Represents a vector with entries of type dai::Real.
 typedef TProb<Real> Prob;
 
 
index 433b7fc..1cbfd06 100644 (file)
@@ -11,8 +11,7 @@
 
 
 /// \file
-/// \brief Defines the SmallSet<T> class, which represents a set; the implementation is optimized for a small number of elements.
-/// \todo Check documentation
+/// \brief Defines the SmallSet<<em>T</em>> class, which represents a set (optimized for a small number of elements).
 
 
 #ifndef __defined_libdai_smallset_h
@@ -27,8 +26,8 @@ namespace dai {
 
 
 /// Represents a set; the implementation is optimized for a small number of elements.
-/** SmallSet uses an ordered std::vector<T> to represent a set; this is faster than
- *  using a std::set<T> if the number of elements is small.
+/** SmallSet uses an ordered <tt>std::vector<</tt><em>T</em><tt>></tt> to represent a set; this is faster than
+ *  using a <tt>std::set<</tt><em>T</em><tt>></tt> if the number of elements is small.
  *  \tparam T Should be less-than-comparable.
  */
 template <typename T>
@@ -40,15 +39,15 @@ class SmallSet {
     public:
     /// @name Constructors and destructors
     //@{
-        /// Default constructor (construct an empty set)
+        /// Default constructor (constructs an empty set)
         SmallSet() : _elements() {}
 
-        /// Construct a set with one element
+        /// Construct a set consisting of one element
         SmallSet( const T &t ) : _elements() {
             _elements.push_back( t );
         }
 
-        /// Construct a set with two elements
+        /// Construct a set consisting of two elements
         SmallSet( const T &t1, const T &t2 ) {
             if( t1 < t2 ) {
                 _elements.push_back( t1 );
@@ -61,10 +60,10 @@ class SmallSet {
         }
 
         /// Construct a SmallSet from a range of elements.
-        /** \tparam TIterator Iterates over instances of type T.
+        /** \tparam TIterator Iterates over instances of type \a T.
          *  \param begin Points to first element to be added.
          *  \param end Points just beyond last element to be added.
-         *  \param sizeHint For efficiency, the number of elements can be speficied by sizeHint.
+         *  \param sizeHint For efficiency, the number of elements can be speficied by \a sizeHint.
          */
         template <typename TIterator>
         SmallSet( TIterator begin, TIterator end, size_t sizeHint=0 ) {
@@ -78,28 +77,28 @@ class SmallSet {
 
     /// @name Operators for set-theoretic operations
     //@{
-        /// Set-minus operator: returns all elements in *this, except those in x
+        /// Set-minus operator: returns all elements in \c *this, except those in \a x
         SmallSet operator/ ( const SmallSet& x ) const {
             SmallSet res;
             std::set_difference( _elements.begin(), _elements.end(), x._elements.begin(), x._elements.end(), inserter( res._elements, res._elements.begin() ) );
             return res;
         }
 
-        /// Set-union operator: returns all elements in *this, plus those in x
+        /// Set-union operator: returns all elements in \c *this, plus those in \a x
         SmallSet operator| ( const SmallSet& x ) const {
             SmallSet res;
             std::set_union( _elements.begin(), _elements.end(), x._elements.begin(), x._elements.end(), inserter( res._elements, res._elements.begin() ) );
             return res;
         }
 
-        /// Set-intersection operator: returns all elements in *this that are also contained in x
+        /// Set-intersection operator: returns all elements in \c *this that are also contained in \a x
         SmallSet operator& ( const SmallSet& x ) const {
             SmallSet res;
             std::set_intersection( _elements.begin(), _elements.end(), x._elements.begin(), x._elements.end(), inserter( res._elements, res._elements.begin() ) );
             return res;
         }
 
-        /// Erases from *this all elements in x
+        /// Erases from \c *this all elements in \a x
         SmallSet& operator/= ( const SmallSet& x ) {
             return (*this = (*this / x));
         }
@@ -113,7 +112,7 @@ class SmallSet {
             return *this;
         }
 
-        /// Adds to *this all elements in x
+        /// Adds to \c *this all elements in \a x
         SmallSet& operator|= ( const SmallSet& x ) {
             return( *this = (*this | x) );
         }
@@ -126,17 +125,17 @@ class SmallSet {
             return *this;
         }
 
-        /// Erases from *this all elements not in x
+        /// Erases from \c *this all elements not in \a x
         SmallSet& operator&= ( const SmallSet& x ) {
             return (*this = (*this & x));
         }
 
-        /// Returns true if *this is a subset of x
+        /// Returns \c true if \c *this is a subset of \a x
         bool operator<< ( const SmallSet& x ) const {
             return std::includes( x._elements.begin(), x._elements.end(), _elements.begin(), _elements.end() );
         }
 
-        /// Returns true if x is a subset of *this
+        /// Returns \c true if \a x is a subset of \c *this
         bool operator>> ( const SmallSet& x ) const {
             return std::includes( _elements.begin(), _elements.end(), x._elements.begin(), x._elements.end() );
         }
@@ -144,12 +143,12 @@ class SmallSet {
 
     /// @name Queries
     //@{
-        /// Returns true if *this and x have elements in common
+        /// Returns \c true if \c *this and \a x have elements in common
         bool intersects( const SmallSet& x ) const {
             return( (*this & x).size() > 0 );
         }
 
-        /// Returns true if *this contains the element t
+        /// Returns \c true if \c *this contains the element \a t
         bool contains( const T &t ) const {
             return std::binary_search( _elements.begin(), _elements.end(), t );
         }
@@ -157,7 +156,7 @@ class SmallSet {
         /// Returns number of elements
         typename std::vector<T>::size_type size() const { return _elements.size(); }
 
-        /// Returns whether the SmallSet is empty
+        /// Returns whether \c *this is empty
         bool empty() const { return _elements.size() == 0; }
     //@}
 
@@ -195,12 +194,12 @@ class SmallSet {
 
     /// @name Comparison operators
     //@{
-        /// Returns true if a and b are identical
+        /// Returns \c true if \a a and \a b are identical
         friend bool operator==( const SmallSet &a, const SmallSet &b ) {
             return (a._elements == b._elements);
         }
 
-        /// Returns true if a and b are not identical
+        /// Returns \c true if \a a and \a b are not identical
         friend bool operator!=( const SmallSet &a, const SmallSet &b ) {
             return !(a._elements == b._elements);
         }
index 1952562..480c3ef 100644 (file)
 #define DAI_IFVERB(n, stmt) if(props.verbose>=n) { cerr << stmt; }
 
 
-/// Real number (alias for double, which could be changed to long double if necessary)
-typedef double Real;
-
-
 #ifdef WINDOWS
     /// Returns true if argument is NAN (Not A Number)
     bool isnan( double x );
@@ -90,6 +86,20 @@ typedef double Real;
 namespace dai {
 
 
+/// Real number (alias for double, which could be changed to long double if necessary)
+typedef double Real;
+
+/// Returns logarithm of \a x
+inline Real log( Real x ) {
+    return std::log(x);
+}
+
+/// Returns exponent of \a x
+inline Real exp( Real x ) {
+    return std::exp(x);
+}
+
+
 #ifdef WINDOWS
     /// hash_map is an alias for std::map.
     /** Since there is no TR1 unordered_map implementation available yet, we fall back on std::map.
@@ -109,14 +119,21 @@ namespace dai {
 double toc();
 
 
+/// Returns absolute value of \a t
+template<class T>
+inline T abs( const T &t ) {
+    return (t < 0) ? (-t) : t;
+}
+
+
 /// Sets the random seed
 void rnd_seed( size_t seed );
 
 /// Returns a real number, distributed uniformly on [0,1)
-double rnd_uniform();
+Real rnd_uniform();
 
 /// Returns a real number from a standard-normal distribution
-double rnd_stdnormal();
+Real rnd_stdnormal();
 
 /// Returns a random integer in interval [min, max]
 int rnd_int( int min, int max );
index 7232615..87fc45b 100644 (file)
@@ -24,6 +24,26 @@ using namespace std;
 const char *HAK::Name = "HAK";
 
 
+
+/// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
+template <class T>
+TFactor<T>& makePositive( TFactor<T> &f, T epsilon ) {
+    for( size_t t = 0; t < f.states(); t++ )
+        if( (0 < f[t]) && (f[t] < epsilon) )
+            f[t] = epsilon;
+    return f;
+}
+
+/// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
+template <class T>
+TFactor<T>& makeZero( TFactor<T> &f, T epsilon ) {
+    for( size_t t = 0; t < f.states(); t++ )
+        if( f[t] < epsilon && f[t] > -epsilon )
+            f[t] = 0;
+    return f;
+}
+
+
 void HAK::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("tol") );
     DAI_ASSERT( opts.hasKey("maxiter") );
index 3b2f0a0..73fdf8f 100644 (file)
@@ -69,27 +69,27 @@ typedef boost::minstd_rand _rnd_gen_type;  // Try boost::mt19937 or boost::ecuye
 _rnd_gen_type _rnd_gen(42U);
 
 /// Uniform distribution with values between 0 and 1 (0 inclusive, 1 exclusive).
-boost::uniform_real<> _uni_dist(0,1);
+boost::uniform_real<Real> _uni_dist(0,1);
 
 /// Global uniform random random number
-boost::variate_generator<_rnd_gen_type&, boost::uniform_real<> > _uni_rnd(_rnd_gen, _uni_dist);
+boost::variate_generator<_rnd_gen_type&, boost::uniform_real<Real> > _uni_rnd(_rnd_gen, _uni_dist);
 
 /// Normal distribution with mean 0 and standard deviation 1.
 boost::normal_distribution<> _normal_dist;
 
 /// Global random number generator with standard normal distribution
-boost::variate_generator<_rnd_gen_type&, boost::normal_distribution<> > _normal_rnd(_rnd_gen, _normal_dist);
+boost::variate_generator<_rnd_gen_type&, boost::normal_distribution<Real> > _normal_rnd(_rnd_gen, _normal_dist);
 
 
 void rnd_seed( size_t seed ) {
     _rnd_gen.seed(seed);
 }
 
-double rnd_uniform() {
+Real rnd_uniform() {
     return _uni_rnd();
 }
 
-double rnd_stdnormal() {
+Real rnd_stdnormal() {
     return _normal_rnd();
 }
 
index 516908a..f495424 100644 (file)
@@ -44,8 +44,8 @@ typedef matrix::index_array_type::const_iterator matrix_icit;
 Factor BinaryFactor( const Var &n, double field ) {
     DAI_ASSERT( n.states() == 2 );
     double buf[2];
-    buf[0] = exp(-field);
-    buf[1] = exp(field);
+    buf[0] = std::exp(-field);
+    buf[1] = std::exp(field);
     return Factor(n, &buf[0]);
 }
 
@@ -55,8 +55,8 @@ Factor BinaryFactor( const Var &n1, const Var &n2, double coupling ) {
     DAI_ASSERT( n2.states() == 2 );
     DAI_ASSERT( n1 != n2 );
     double buf[4];
-    buf[0] = (buf[3] = exp(coupling));
-    buf[1] = (buf[2] = exp(-coupling));
+    buf[0] = (buf[3] = std::exp(coupling));
+    buf[1] = (buf[2] = std::exp(-coupling));
     return Factor( VarSet(n1, n2), &buf[0] );
 }
 
@@ -64,7 +64,7 @@ Factor BinaryFactor( const Var &n1, const Var &n2, double coupling ) {
 Factor RandomFactor( const VarSet &ns, double beta ) {
     Factor fac( ns );
     for( size_t t = 0; t < fac.states(); t++ )
-        fac[t] = exp(rnd_stdnormal() * beta);
+        fac[t] = std::exp(rnd_stdnormal() * beta);
     return fac;
 }
 
@@ -73,7 +73,7 @@ Factor PottsFactor( const Var &n1, const Var &n2, double beta ) {
     Factor fac( VarSet( n1, n2 ), 1.0 );
     DAI_ASSERT( n1.states() == n2.states() );
     for( size_t s = 0; s < n1.states(); s++ )
-        fac[ s * (n1.states() + 1) ] = exp(beta);
+        fac[ s * (n1.states() + 1) ] = std::exp(beta);
     return fac;
 }