Improved prob.h/cpp code:
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 6 Apr 2010 08:53:14 +0000 (10:53 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 6 Apr 2010 08:53:14 +0000 (10:53 +0200)
- Deprecated TProb::TProb( TIterator begin, TIterator end, size_t sizeHint=0 )
  constructor: the sizeHint argument no longer defaults to 0
- Deprecated TProb::accumulate( T init, binOp op1, unOp op2 );
  instead, TProb::accumulateSum( T init, unOp op ) and
  TProb::accumulateMax( T init, unOp op, bool minimize ) should be used.
- Deprecated TProb::NormType and TProb::DistType;
  ProbNormType and ProbDistType should be used instead.

15 files changed:
ChangeLog
include/dai/factor.h
include/dai/prob.h
include/dai/util.h
src/bbp.cpp
src/bp.cpp
src/cbp.cpp
src/fbp.cpp
src/hak.cpp
src/lc.cpp
src/mf.cpp
src/treeep.cpp
tests/testdai.cpp
tests/unit/factor.cpp
tests/unit/prob.cpp

index 8d185e8..9a81f7e 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -20,6 +20,13 @@ TODO:
     equivalent to "operator[](size_t) const" and set() should
     be used instead of the non-const operator[], which has been deprecated
 * Improved prob.h/cpp:
+  - Deprecated TProb::TProb( TIterator begin, TIterator end, size_t sizeHint=0 )
+    constructor: the sizeHint argument no longer defaults to 0
+  - Deprecated TProb::accumulate( T init, binOp op1, unOp op2 );
+    instead, TProb::accumulateSum( T init, unOp op ) and
+    TProb::accumulateMax( T init, unOp op, bool minimize ) should be used.
+  - Deprecated TProb::NormType and TProb::DistType;
+    ProbNormType and ProbDistType should be used instead.
   - Added TProb<T>::operator==( const TProb<T> &q )
   - TProb<T>::accumulate() now also applies op2 to init
   - Fixed bug by renaming TProb<T>::operator<=() to TProb<T>::operator<()
@@ -28,7 +35,7 @@ TODO:
   - Added get(size_t) and set(size_t,T) operators; get() is 
     equivalent to "operator[](size_t) const" and set() should
     be used instead of the non-const operator[], which has been deprecated
-  - Added TProb<T>::resize(size_t, T)
+  - Added TProb<T>::resize(size_t)
 * Improved index.h/cpp:
   - Added multifor::reset()
 * Improved properties.h/cpp:
index 9925835..ca596bb 100644 (file)
@@ -55,12 +55,13 @@ namespace dai {
  *  \todo Define a better fileformat for .fg files (maybe using XML)?
  *  \todo Add support for sparse factors.
  */
-template <typename T> class TFactor {
+template <typename T>
+class TFactor {
     private:
         /// Stores the variables on which the factor depends
-        VarSet      _vs;
+        VarSet _vs;
         /// Stores the factor values
-        TProb<T>    _p;
+        TProb<T> _p;
 
     public:
     /// \name Constructors and destructors
@@ -110,6 +111,15 @@ template <typename T> class TFactor {
         }
     //@}
 
+    /// \name Get/set individual entries
+    //@{
+        /// Sets \a i 'th entry to \a val
+        void set( size_t i, T val ) { _p.set( i, val ); }
+
+        /// Gets \a i 'th entry
+        T get( size_t i ) const { return _p[i]; }
+    //@}
+
     /// \name Queries
     //@{
         /// Returns constant reference to value vector
@@ -125,12 +135,6 @@ template <typename T> class TFactor {
         /// \deprecated Please use dai::TFactor::set() instead
         T& operator[] (size_t i) { return _p[i]; }
 
-        /// Gets \a i 'th entry of the value vector
-        T get( size_t i ) const { return _p[i]; }
-
-        /// Sets \a i 'th entry of the value vector to \a val
-        void set( size_t i, T val ) { _p.set( i, val ); }
-
         /// Returns constant reference to variable set (i.e., the variables on which the factor depends)
         const VarSet& vars() const { return _vs; }
 
@@ -219,7 +223,7 @@ template <typename T> class TFactor {
             x._p = _p.log(zero);
             return x;
         }
-        
+
         /// Returns pointwise inverse
         /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
          */
@@ -233,7 +237,7 @@ template <typename T> class TFactor {
         /// Returns normalized copy of \c *this, using the specified norm
         /** \throw NOT_NORMALIZABLE if the norm is zero
          */
-        TFactor<T> normalized( typename TProb<T>::NormType norm=TProb<T>::NORMPROB ) const {
+        TFactor<T> normalized( ProbNormType norm=NORMPROB ) const {
             TFactor<T> x;
             x._vs = _vs;
             x._p = _p.normalized( norm );
@@ -263,7 +267,7 @@ template <typename T> class TFactor {
         /// Normalizes factor using the specified norm
         /** \throw NOT_NORMALIZABLE if the norm is zero
          */
-        T normalize( typename TProb<T>::NormType norm=TProb<T>::NORMPROB ) { return _p.normalize( norm ); }
+        T normalize( ProbNormType norm=NORMPROB ) { return _p.normalize( norm ); }
     //@}
 
     /// \name Operations with scalars
@@ -413,7 +417,7 @@ template <typename T> class TFactor {
                 result._p.p().clear();
                 result._p.p().reserve( N );
                 for( size_t i = 0; i < N; i++, ++i_f, ++i_g )
-                    result._p.p().push_back( op( _p[i_f], g._p[i_g] ) );
+                    result._p.p().push_back( op( _p[i_f], g[i_g] ) );
             }
             return result;
         }
@@ -520,7 +524,7 @@ template<typename T> TFactor<T> TFactor<T>::marginal(const VarSet &vars, bool no
         res.set( i_res, res[i_res] + _p[i] );
 
     if( normed )
-        res.normalize( TProb<T>::NORMPROB );
+        res.normalize( NORMPROB );
 
     return res;
 }
@@ -537,7 +541,7 @@ template<typename T> TFactor<T> TFactor<T>::maxMarginal(const VarSet &vars, bool
             res.set( i_res, _p[i] );
 
     if( normed )
-        res.normalize( TProb<T>::NORMPROB );
+        res.normalize( NORMPROB );
 
     return res;
 }
@@ -588,7 +592,7 @@ template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T
 /** \relates TFactor
  *  \pre f.vars() == g.vars()
  */
-template<typename T> T dist( const TFactor<T> &f, const TFactor<T> &g, typename TProb<T>::DistType dt ) {
+template<typename T> T dist( const TFactor<T> &f, const TFactor<T> &g, ProbDistType dt ) {
     if( f.vars().empty() || g.vars().empty() )
         return -1;
     else {
@@ -627,7 +631,7 @@ template<typename T> T MutualInfo(const TFactor<T> &f) {
     VarSet::const_iterator it = f.vars().begin();
     Var i = *it; it++; Var j = *it;
     TFactor<T> projection = f.marginal(i) * f.marginal(j);
-    return dist( f.normalized(), projection, TProb<T>::DISTKL );
+    return dist( f.normalized(), projection, DISTKL );
 }
 
 
index a44da76..7212f7a 100644 (file)
@@ -198,6 +198,7 @@ class TProb {
     public:
         /// Type of data structure used for storing the values
         typedef std::vector<T> container_type;
+        typedef TProb<T> this_type;
 
     private:
         /// The data structure that stores the values
@@ -208,6 +209,7 @@ class TProb {
         /**
          *  - NORMPROB means that the sum of all entries should be 1;
          *  - NORMLINF means that the maximum absolute value of all entries should be 1.
+         *  \deprecated Please use dai::ProbNormType instead.
          */
         typedef enum { NORMPROB, NORMLINF } NormType;
         /// Enumerates different distance measures between probability measures.
@@ -217,6 +219,7 @@ class TProb {
          *  - 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$).
          *  - DISTHEL is the Hellinger distance (\f$\frac{1}{2}\sum_i (\sqrt{p_i}-\sqrt{q_i})^2\f$).
+         *  \deprecated Please use dai::ProbDistType instead.
          */
         typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL, DISTHEL } DistType;
 
@@ -235,7 +238,9 @@ class TProb {
         /** \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 \a sizeHint.
+         *  \param sizeHint For efficiency, the number of entries can be speficied by \a sizeHint;
+         *    the value 0 can be given if the size is unknown, but this will result in a performance penalty.
+         *  \deprecated In future libDAI versions, the \a sizeHint argument will no longer default to 0.
          */
         template <typename TIterator>
         TProb( TIterator begin, TIterator end, size_t sizeHint=0 ) : _p() {
@@ -288,12 +293,12 @@ class TProb {
 
     /// \name Miscellaneous operations
     //@{
-        void resize( size_t sz, T c = T() ) {
-            _p.resize( sz, c );
+        void resize( size_t sz ) {
+            _p.resize( sz );
         }
     //@}
 
-    /// \name Queries
+    /// \name Get/set individual entries
     //@{
         /// Gets \a i 'th entry
         T get( size_t i ) const { 
@@ -309,11 +314,14 @@ class TProb {
             DAI_DEBASSERT( i < _p.size() );
             _p[i] = val;
         }
+    //@}
 
-        /// Returns a const reference to the wrapped vector
+    /// \name Queries
+    //@{
+        /// Returns a const reference to the wrapped container
         const container_type& p() const { return _p; }
 
-        /// Returns a reference to the wrapped vector
+        /// Returns a reference to the wrapped container
         container_type& p() { return _p; }
 
         /// Returns a copy of the \a i 'th entry
@@ -335,6 +343,7 @@ class TProb {
          *      t = op1( t, op2(*it) );
          *  return t;
          *  \endcode
+         *  \deprecated Please use dai::TProb::accumulateSum or dai::TProb::accumulateMax instead
          */
         template<typename binOp, typename unOp> T accumulate( T init, binOp op1, unOp op2 ) const {
             T t = op2(init);
@@ -343,23 +352,61 @@ class TProb {
             return t;
         }
 
+
+        /// Accumulate all values (similar to std::accumulate) by summing
+        /** The following calculation is done:
+         *  \code
+         *  T t = op(init);
+         *  for( const_iterator it = begin(); it != end(); it++ )
+         *      t += op(*it);
+         *  return t;
+         *  \endcode
+         */
+        template<typename unOp> T accumulateSum( T init, unOp op ) const {
+            T t = op(init);
+            for( const_iterator it = begin(); it != end(); it++ )
+                t += op(*it);
+            return t;
+        }
+
+        /// Accumulate all values (similar to std::accumulate) by maximization/minimization
+        /** The following calculation is done (with "max" replaced by "min" if \a minimize == \c true):
+         *  \code
+         *  T t = op(init);
+         *  for( const_iterator it = begin(); it != end(); it++ )
+         *      t = std::max( t, op(*it) );
+         *  return t;
+         *  \endcode
+         */
+        template<typename unOp> T accumulateMax( T init, unOp op, bool minimize ) const {
+            T t = op(init);
+            if( minimize ) {
+                for( const_iterator it = begin(); it != end(); it++ )
+                    t = std::min( t, op(*it) );
+            } else {
+                for( const_iterator it = begin(); it != end(); it++ )
+                    t = std::max( t, op(*it) );
+            }
+            return t;
+        }
+
         /// Returns the Shannon entropy of \c *this, \f$-\sum_i p_i \log p_i\f$
-        T entropy() const { return -accumulate( (T)0, std::plus<T>(), fo_plog0p<T>() ); }
+        T entropy() const { return -accumulateSum( (T)0, fo_plog0p<T>() ); }
 
         /// Returns maximum value of all entries
-        T max() const { return accumulate( (T)(-INFINITY), fo_max<T>(), fo_id<T>() ); }
+        T max() const { return accumulateMax( (T)(-INFINITY), fo_id<T>(), false ); }
 
         /// Returns minimum value of all entries
-        T min() const { return accumulate( (T)INFINITY, fo_min<T>(), fo_id<T>() ); }
+        T min() const { return accumulateMax( (T)INFINITY, fo_id<T>(), true ); }
 
         /// Returns sum of all entries
-        T sum() const { return accumulate( (T)0, std::plus<T>(), fo_id<T>() ); }
+        T sum() const { return accumulateSum( (T)0, fo_id<T>() ); }
 
         /// Return sum of absolute value of all entries
-        T sumAbs() const { return accumulate( (T)0, std::plus<T>(), fo_abs<T>() ); }
+        T sumAbs() const { return accumulateSum( (T)0, fo_abs<T>() ); }
 
         /// Returns maximum absolute value of all entries
-        T maxAbs() const { return accumulate( (T)0, fo_max<T>(), fo_abs<T>() ); }
+        T maxAbs() const { return accumulateMax( (T)0, fo_abs<T>(), false ); }
 
         /// Returns \c true if one or more entries are NaN
         bool hasNaNs() const {
@@ -405,13 +452,13 @@ class TProb {
         /// Lexicographical comparison
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        bool operator<( const TProb<T>& q ) const {
+        bool operator<( const this_type& q ) const {
             DAI_DEBASSERT( size() == q.size() );
             return lexicographical_compare( begin(), end(), q.begin(), q.end() );
         }
 
         /// Comparison
-        bool operator==( const TProb<T>& q ) const {
+        bool operator==( const this_type& q ) const {
             if( size() != q.size() )
                 return false;
             return p() == q.p();
@@ -421,26 +468,26 @@ class TProb {
     /// \name Unary transformations
     //@{
         /// Returns the result of applying operation \a op pointwise on \c *this
-        template<typename unaryOp> TProb<T> pwUnaryTr( unaryOp op ) const {
-            TProb<T> r;
+        template<typename unaryOp> this_type pwUnaryTr( unaryOp op ) const {
+            this_type r;
             r._p.reserve( size() );
             std::transform( _p.begin(), _p.end(), back_inserter( r._p ), op );
             return r;
         }
 
         /// Returns negative of \c *this
-        TProb<T> operator- () const { return pwUnaryTr( std::negate<T>() ); }
+        this_type operator- () const { return pwUnaryTr( std::negate<T>() ); }
 
         /// Returns pointwise absolute value
-        TProb<T> abs() const { return pwUnaryTr( fo_abs<T>() ); }
+        this_type abs() const { return pwUnaryTr( fo_abs<T>() ); }
 
         /// Returns pointwise exponent
-        TProb<T> exp() const { return pwUnaryTr( fo_exp<T>() ); }
+        this_type exp() const { return pwUnaryTr( fo_exp<T>() ); }
 
         /// Returns pointwise logarithm
         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
          */
-        TProb<T> log(bool zero=false) const {
+        this_type log(bool zero=false) const {
             if( zero )
                 return pwUnaryTr( fo_log0<T>() );
             else
@@ -450,7 +497,7 @@ class TProb {
         /// Returns pointwise inverse
         /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
          */
-        TProb<T> inverse(bool zero=true) const {
+        this_type inverse(bool zero=true) const {
             if( zero )
                 return pwUnaryTr( fo_inv0<T>() );
             else
@@ -460,11 +507,11 @@ class TProb {
         /// Returns normalized copy of \c *this, using the specified norm
         /** \throw NOT_NORMALIZABLE if the norm is zero
          */
-        TProb<T> normalized( NormType norm = NORMPROB ) const {
+        this_type normalized( ProbNormType norm = dai::NORMPROB ) const {
             T Z = 0;
-            if( norm == NORMPROB )
+            if( norm == dai::NORMPROB )
                 Z = sum();
-            else if( norm == NORMLINF )
+            else if( norm == dai::NORMLINF )
                 Z = maxAbs();
             if( Z == (T)0 ) {
                 DAI_THROW(NOT_NORMALIZABLE);
@@ -477,33 +524,33 @@ class TProb {
     /// \name Unary operations
     //@{
         /// Applies unary operation \a op pointwise
-        template<typename unaryOp> TProb<T>& pwUnaryOp( unaryOp op ) {
+        template<typename unaryOp> this_type& pwUnaryOp( unaryOp op ) {
             std::transform( _p.begin(), _p.end(), _p.begin(), op );
             return *this;
         }
 
         /// Draws all entries i.i.d. from a uniform distribution on [0,1)
-        TProb<T>& randomize() {
+        this_type& randomize() {
             std::generate( _p.begin(), _p.end(), rnd_uniform );
             return *this;
         }
 
         /// Sets all entries to \f$1/n\f$ where \a n is the length of the vector
-        TProb<T>& setUniform () {
+        this_type& setUniform () {
             fill( (T)1 / size() );
             return *this;
         }
 
         /// Applies absolute value pointwise
-        TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
+        this_type& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
 
         /// Applies exponent pointwise
-        TProb<T>& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
+        this_type& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
 
         /// Applies logarithm pointwise
         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
          */
-        TProb<T>& takeLog(bool zero=false) {
+        this_type& takeLog(bool zero=false) {
             if( zero ) {
                 return pwUnaryOp( fo_log0<T>() );
             } else
@@ -513,11 +560,11 @@ class TProb {
         /// Normalizes vector using the specified norm
         /** \throw NOT_NORMALIZABLE if the norm is zero
          */
-        T normalize( NormType norm=NORMPROB ) {
+        T normalize( ProbNormType norm=dai::NORMPROB ) {
             T Z = 0;
-            if( norm == NORMPROB )
+            if( norm == dai::NORMPROB )
                 Z = sum();
-            else if( norm == NORMLINF )
+            else if( norm == dai::NORMLINF )
                 Z = maxAbs();
             if( Z == (T)0 )
                 DAI_THROW(NOT_NORMALIZABLE);
@@ -530,13 +577,13 @@ class TProb {
     /// \name Operations with scalars
     //@{
         /// Sets all entries to \a x
-        TProb<T> & fill( T x ) {
+        this_type& fill( T x ) {
             std::fill( _p.begin(), _p.end(), x );
             return *this;
         }
 
         /// Adds scalar \a x to each entry
-        TProb<T>& operator+= (T x) {
+        this_type& operator+= (T x) {
             if( x != 0 )
                 return pwUnaryOp( std::bind2nd( std::plus<T>(), x ) );
             else
@@ -544,7 +591,7 @@ class TProb {
         }
 
         /// Subtracts scalar \a x from each entry
-        TProb<T>& operator-= (T x) {
+        this_type& operator-= (T x) {
             if( x != 0 )
                 return pwUnaryOp( std::bind2nd( std::minus<T>(), x ) );
             else
@@ -552,7 +599,7 @@ class TProb {
         }
 
         /// Multiplies each entry with scalar \a x
-        TProb<T>& operator*= (T x) {
+        this_type& operator*= (T x) {
             if( x != 1 )
                 return pwUnaryOp( std::bind2nd( std::multiplies<T>(), x ) );
             else
@@ -560,7 +607,7 @@ class TProb {
         }
 
         /// Divides each entry by scalar \a x, where division by 0 yields 0
-        TProb<T>& operator/= (T x) {
+        this_type& operator/= (T x) {
             if( x != 1 )
                 return pwUnaryOp( std::bind2nd( fo_divides0<T>(), x ) );
             else
@@ -568,7 +615,7 @@ class TProb {
         }
 
         /// Raises entries to the power \a x
-        TProb<T>& operator^= (T x) {
+        this_type& operator^= (T x) {
             if( x != (T)1 )
                 return pwUnaryOp( std::bind2nd( fo_pow<T>(), x) );
             else
@@ -579,19 +626,19 @@ class TProb {
     /// \name Transformations with scalars
     //@{
         /// Returns sum of \c *this and scalar \a x
-        TProb<T> operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }
+        this_type operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }
 
         /// Returns difference of \c *this and scalar \a x
-        TProb<T> operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }
+        this_type operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }
 
         /// Returns product of \c *this with scalar \a x
-        TProb<T> operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }
+        this_type operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }
 
         /// Returns quotient of \c *this and scalar \a x, where division by 0 yields 0
-        TProb<T> operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }
+        this_type operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }
 
         /// Returns \c *this raised to the power \a x
-        TProb<T> operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }
+        this_type operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }
     //@}
 
     /// \name Operations with other equally-sized vectors
@@ -601,7 +648,7 @@ class TProb {
          *  \param q Right operand
          *  \param op Operation of type \a binaryOp
          */
-        template<typename binaryOp> TProb<T>& pwBinaryOp( const TProb<T> &q, binaryOp op ) {
+        template<typename binaryOp> this_type& pwBinaryOp( const this_type &q, binaryOp op ) {
             DAI_DEBASSERT( size() == q.size() );
             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), op );
             return *this;
@@ -610,34 +657,34 @@ class TProb {
         /// Pointwise addition with \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T>& operator+= (const TProb<T> & q) { return pwBinaryOp( q, std::plus<T>() ); }
+        this_type& operator+= (const this_type & q) { return pwBinaryOp( q, std::plus<T>() ); }
 
         /// Pointwise subtraction of \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T>& operator-= (const TProb<T> & q) { return pwBinaryOp( q, std::minus<T>() ); }
+        this_type& operator-= (const this_type & q) { return pwBinaryOp( q, std::minus<T>() ); }
 
         /// Pointwise multiplication with \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T>& operator*= (const TProb<T> & q) { return pwBinaryOp( q, std::multiplies<T>() ); }
+        this_type& operator*= (const this_type & q) { return pwBinaryOp( q, std::multiplies<T>() ); }
 
         /// Pointwise division by \a q, where division by 0 yields 0
         /** \pre <tt>this->size() == q.size()</tt>
          *  \see divide(const TProb<T> &)
          */
-        TProb<T>& operator/= (const TProb<T> & q) { return pwBinaryOp( q, fo_divides0<T>() ); }
+        this_type& operator/= (const this_type & q) { return pwBinaryOp( q, fo_divides0<T>() ); }
 
         /// Pointwise division by \a q, where division by 0 yields +Inf
         /** \pre <tt>this->size() == q.size()</tt>
          *  \see operator/=(const TProb<T> &)
          */
-        TProb<T>& divide (const TProb<T> & q) { return pwBinaryOp( q, std::divides<T>() ); }
+        this_type& divide (const this_type & q) { return pwBinaryOp( q, std::divides<T>() ); }
 
         /// Pointwise power
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T>& operator^= (const TProb<T> & q) { return pwBinaryOp( q, fo_pow<T>() ); }
+        this_type& operator^= (const this_type & q) { return pwBinaryOp( q, fo_pow<T>() ); }
     //@}
 
     /// \name Transformations with other equally-sized vectors
@@ -647,7 +694,7 @@ class TProb {
          *  \param q Right operand
          *  \param op Operation of type \a binaryOp
          */
-        template<typename binaryOp> TProb<T> pwBinaryTr( const TProb<T> &q, binaryOp op ) const {
+        template<typename binaryOp> this_type pwBinaryTr( const this_type &q, binaryOp op ) const {
             DAI_DEBASSERT( size() == q.size() );
             TProb<T> r;
             r._p.reserve( size() );
@@ -658,40 +705,40 @@ class TProb {
         /// Returns sum of \c *this and \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T> operator+ ( const TProb<T>& q ) const { return pwBinaryTr( q, std::plus<T>() ); }
+        this_type operator+ ( const this_type& q ) const { return pwBinaryTr( q, std::plus<T>() ); }
 
         /// Return \c *this minus \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T> operator- ( const TProb<T>& q ) const { return pwBinaryTr( q, std::minus<T>() ); }
+        this_type operator- ( const this_type& q ) const { return pwBinaryTr( q, std::minus<T>() ); }
 
         /// Return product of \c *this with \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T> operator* ( const TProb<T> &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }
+        this_type operator* ( const this_type &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }
 
         /// Returns quotient of \c *this with \a q, where division by 0 yields 0
         /** \pre <tt>this->size() == q.size()</tt>
          *  \see divided_by(const TProb<T> &)
          */
-        TProb<T> operator/ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }
+        this_type operator/ ( const this_type &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }
 
         /// Pointwise division by \a q, where division by 0 yields +Inf
         /** \pre <tt>this->size() == q.size()</tt>
          *  \see operator/(const TProb<T> &)
          */
-        TProb<T> divided_by( const TProb<T> &q ) const { return pwBinaryTr( q, std::divides<T>() ); }
+        this_type divided_by( const this_type &q ) const { return pwBinaryTr( q, std::divides<T>() ); }
 
         /// Returns \c *this to the power \a q
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        TProb<T> operator^ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }
+        this_type operator^ ( const this_type &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }
     //@}
 
         /// Performs a generalized inner product, similar to std::inner_product
         /** \pre <tt>this->size() == q.size()</tt>
          */
-        template<typename binOp1, typename binOp2> T innerProduct( const TProb<T> &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
+        template<typename binOp1, typename binOp2> T innerProduct( const this_type &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
             DAI_DEBASSERT( size() == q.size() );
             return std::inner_product( begin(), end(), q.begin(), init, binaryOp1, binaryOp2 );
         }
@@ -702,17 +749,17 @@ class TProb {
 /** \relates TProb
  *  \pre <tt>this->size() == q.size()</tt>
  */
-template<typename T> T 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, ProbDistType dt ) {
     switch( dt ) {
-        case TProb<T>::DISTL1:
+        case DISTL1:
             return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() );
-        case TProb<T>::DISTLINF:
+        case DISTLINF:
             return p.innerProduct( q, (T)0, fo_max<T>(), fo_absdiff<T>() );
-        case TProb<T>::DISTTV:
+        case DISTTV:
             return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() ) / 2;
-        case TProb<T>::DISTKL:
+        case DISTKL:
             return p.innerProduct( q, (T)0, std::plus<T>(), fo_KL<T>() );
-        case TProb<T>::DISTHEL:
+        case DISTHEL:
             return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
         default:
             DAI_THROW(UNKNOWN_ENUM_VALUE);
@@ -726,8 +773,8 @@ template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TPro
  */
 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
     os << "(";
-    for( typename TProb<T>::const_iterator it = p.begin(); it != p.end(); it++ )
-        os << (it != p.begin() ? ", " : "") << *it;
+    for( size_t i = 0; i < p.size(); i++ )
+        os << ((i != 0) ? ", " : "") << p.get(i);
     os << ")";
     return os;
 }
index 8fe4f07..a5bc03c 100644 (file)
@@ -186,6 +186,24 @@ std::vector<T> concat( const std::vector<T>& u, const std::vector<T>& v ) {
 void tokenizeString( const std::string& s, std::vector<std::string>& outTokens, const std::string& delim="\t\n" );
 
 
+/// Enumerates different ways of normalizing a probability measure.
+/**
+ *  - NORMPROB means that the sum of all entries should be 1;
+ *  - NORMLINF means that the maximum absolute value of all entries should be 1.
+ */
+typedef enum { NORMPROB, NORMLINF } ProbNormType;
+
+/// Enumerates different distance measures between probability measures.
+/**
+ *  - 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$).
+ *  - DISTHEL is the Hellinger distance (\f$\frac{1}{2}\sum_i (\sqrt{p_i}-\sqrt{q_i})^2\f$).
+ */
+typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL, DISTHEL } ProbDistType;
+
+
 } // end of namespace dai
 
 
index 8c73591..5a4049d 100644 (file)
@@ -1036,7 +1036,7 @@ Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const P
             cout << "i: " << i
                  << ", p_adj_est: " << p_adj_est
                  << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl;
-            d += dist( p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1 );
+            d += dist( p_adj_est, bbp.adj_psi_V(i), DISTL1 );
         }
     }
     /*    if(1) {
@@ -1086,14 +1086,14 @@ Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const P
                 cerr << "i: " << i << ", I: " << I
                      << ", adj_n_est: " << p_adj_n_est
                      << ", bbp.adj_n(i,I): " << bbp.adj_n(i,I) << endl;
-                d += dist(p_adj_n_est, bbp.adj_n(i,I), Prob::DISTL1);
+                d += dist(p_adj_n_est, bbp.adj_n(i,I), DISTL1);
 
                 Prob p_adj_m_est( adj_m_est );
                 // compare this numerical estimate to the BBP estimate; sum the distances
                 cerr << "i: " << i << ", I: " << I
                      << ", adj_m_est: " << p_adj_m_est
                      << ", bbp.adj_m(I,i): " << bbp.adj_m(I,i) << endl;
-                d += dist(p_adj_m_est, bbp.adj_m(I,i), Prob::DISTL1);
+                d += dist(p_adj_m_est, bbp.adj_m(I,i), DISTL1);
             }
         }
     }
@@ -1120,7 +1120,7 @@ Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const P
             cerr << "i: " << i
                  << ", adj_b_V_est: " << p_adj_b_V_est
                  << ", bbp.adj_b_V(i): " << bbp.adj_b_V(i) << endl;
-            d += dist(p_adj_b_V_est, bbp.adj_b_V(i), Prob::DISTL1);
+            d += dist(p_adj_b_V_est, bbp.adj_b_V(i), DISTL1);
         }
     }
     */
index 4f4f35c..ab18715 100644 (file)
@@ -232,7 +232,7 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
 
     // Update the residual if necessary
     if( props.updates == Properties::UpdateType::SEQMAX )
-        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
+        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), DISTLINF ) );
 }
 
 
@@ -316,12 +316,12 @@ Real BP::run() {
         maxDiff = -INFINITY;
         for( size_t i = 0; i < nrVars(); ++i ) {
             Factor b( beliefV(i) );
-            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], DISTLINF ) );
             oldBeliefsV[i] = b;
         }
         for( size_t I = 0; I < nrFactors(); ++I ) {
             Factor b( beliefF(I) );
-            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], DISTLINF ) );
             oldBeliefsF[I] = b;
         }
 
@@ -418,7 +418,7 @@ Real BP::logZ() const {
     for( size_t i = 0; i < nrVars(); ++i )
         sum += (1.0 - nbV(i).size()) * beliefV(i).entropy();
     for( size_t I = 0; I < nrFactors(); ++I )
-        sum -= dist( beliefF(I), factor(I), Prob::DISTKL );
+        sum -= dist( beliefF(I), factor(I), DISTKL );
     return sum;
 }
 
@@ -455,7 +455,7 @@ void BP::updateMessage( size_t i, size_t _I ) {
         else
             message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
         if( props.updates == Properties::UpdateType::SEQMAX )
-            updateResidual( i, _I, dist( newMessage(i,_I), message(i,_I), Prob::DISTLINF ) );
+            updateResidual( i, _I, dist( newMessage(i,_I), message(i,_I), DISTLINF ) );
     }
 }
 
index 0fbbb44..f565133 100644 (file)
@@ -69,7 +69,7 @@ Real logSumExp( Real a, Real b ) {
 Real dist( const vector<Factor> &b1, const vector<Factor> &b2, size_t nv ) {
     Real d = 0.0;
     for( size_t k = 0; k < nv; k++ )
-        d += dist( b1[k], b2[k], Prob::DISTLINF );
+        d += dist( b1[k], b2[k], DISTLINF );
     return d;
 }
 
@@ -388,7 +388,7 @@ bool CBP::chooseNextClampVar( InfAlg *bp, vector<size_t> &clamped_vars_list, siz
                 Real cost = 0;
                 if( doL1 )
                     for( size_t j = 0; j < nrVars(); j++ )
-                        cost += dist( bp->beliefV(j), bp1->beliefV(j), Prob::DISTL1 );
+                        cost += dist( bp->beliefV(j), bp1->beliefV(j), DISTL1 );
                 else
                     cost = props.bbp_cfn.evaluate( *bp1, &state );
                 if( cost > max_cost || win_k == -1 ) {
index bf140b6..979dda7 100644 (file)
@@ -153,7 +153,7 @@ void FBP::calcNewMessage( size_t i, size_t _I ) {
 
     // Update the residual if necessary
     if( props.updates == Properties::UpdateType::SEQMAX )
-        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
+        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), DISTLINF ) );
 }
 
 
index 60d072b..db8b7da 100644 (file)
@@ -380,12 +380,12 @@ Real HAK::doGBP() {
         maxDiff = -INFINITY;
         for( size_t i = 0; i < nrVars(); ++i ) {
             Factor b = beliefV(i);
-            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], DISTLINF ) );
             oldBeliefsV[i] = b;
         }
         for( size_t I = 0; I < nrFactors(); ++I ) {
             Factor b = beliefF(I);
-            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], DISTLINF ) );
             oldBeliefsF[I] = b;
         }
 
@@ -469,12 +469,12 @@ Real HAK::doDoubleLoop() {
         maxDiff = -INFINITY;
         for( size_t i = 0; i < nrVars(); ++i ) {
             Factor b = beliefV(i);
-            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], DISTLINF ) );
             oldBeliefsV[i] = b;
         }
         for( size_t I = 0; I < nrFactors(); ++I ) {
             Factor b = beliefF(I);
-            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], DISTLINF ) );
             oldBeliefsF[I] = b;
         }
 
index 5aa3dbf..5b360aa 100644 (file)
@@ -307,7 +307,7 @@ Real LC::run() {
         // compare new beliefs with old ones
         maxDiff = -INFINITY;
         for( size_t i = 0; i < nrVars(); i++ ) {
-            maxDiff = std::max( maxDiff, dist( beliefV(i), oldBeliefsV[i], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( beliefV(i), oldBeliefsV[i], DISTLINF ) );
             oldBeliefsV[i] = beliefV(i);
         }
 
index a143608..535ceed 100644 (file)
@@ -131,7 +131,7 @@ Real MF::run() {
             if( props.damping != 0.0 )
                 nb = (nb^(1.0 - props.damping)) * (_beliefs[i]^props.damping);
 
-            maxDiff = std::max( maxDiff, dist( nb, _beliefs[i], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( nb, _beliefs[i], DISTLINF ) );
             _beliefs[i] = nb;
         }
 
index 9030b24..9c52cb7 100644 (file)
@@ -99,7 +99,7 @@ TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opt
                             if( piet.vars() >> ij ) {
                                 piet = piet.marginal( ij );
                                 Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
-                                wg[UEdge(i,findVar(*j))] = dist( piet, pietf, Prob::DISTKL );
+                                wg[UEdge(i,findVar(*j))] = dist( piet, pietf, DISTKL );
                             } else
                                 wg[UEdge(i,findVar(*j))] = 0;
                         } else {
@@ -196,7 +196,7 @@ Real TreeEP::run() {
         vector<Factor> newBeliefs = beliefs();
         maxDiff = -INFINITY;
         for( size_t t = 0; t < oldBeliefs.size(); t++ )
-            maxDiff = std::max( maxDiff, dist( newBeliefs[t], oldBeliefs[t], Prob::DISTLINF ) );
+            maxDiff = std::max( maxDiff, dist( newBeliefs[t], oldBeliefs[t], DISTLINF ) );
         swap( newBeliefs, oldBeliefs );
 
         if( props.verbose >= 3 )
index f47ef4d..64f936c 100644 (file)
@@ -153,7 +153,7 @@ class TestDAI {
             err.clear();
             err.reserve( varMarginals.size() );
             for( size_t i = 0; i < varMarginals.size(); i++ )
-                err.push_back( dist( varMarginals[i], x.varMarginals[i], Prob::DISTTV ) );
+                err.push_back( dist( varMarginals[i], x.varMarginals[i], DISTTV ) );
         }
 
         /// Calculate total variation distance of variable marginals with respect to those in \a x
@@ -161,7 +161,7 @@ class TestDAI {
             err.clear();
             err.reserve( varMarginals.size() );
             for( size_t i = 0; i < varMarginals.size(); i++ )
-                err.push_back( dist( varMarginals[i], x[i], Prob::DISTTV ) );
+                err.push_back( dist( varMarginals[i], x[i], DISTTV ) );
         }
 
         /// Return maximum error
index 3beaf4f..5df6991 100644 (file)
@@ -249,13 +249,13 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 0.5 );
 
-    y = x.normalized( Prob::NORMPROB );
+    y = x.normalized( NORMPROB );
     BOOST_CHECK_EQUAL( y[0], 0.5 );
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 0.5 );
 
     x.set( 0, -2.0 );
-    y = x.normalized( Prob::NORMLINF );
+    y = x.normalized( NORMLINF );
     BOOST_CHECK_EQUAL( y[0], -1.0 );
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 1.0 );
@@ -304,14 +304,14 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     BOOST_CHECK( x == y );
 
     x = xorg;
-    BOOST_CHECK_EQUAL( x.normalize( Prob::NORMPROB ), 3.0 );
+    BOOST_CHECK_EQUAL( x.normalize( NORMPROB ), 3.0 );
     BOOST_CHECK( x == y );
 
     y.set( 0, 2.0 / 2.0 );
     y.set( 1, 0.0 / 2.0 );
     y.set( 2, 1.0 / 2.0 );
     x = xorg;
-    BOOST_CHECK_EQUAL( x.normalize( Prob::NORMLINF ), 2.0 );
+    BOOST_CHECK_EQUAL( x.normalize( NORMLINF ), 2.0 );
     BOOST_CHECK( x == y );
 
     xorg.set( 0, -2.0 );
@@ -814,30 +814,30 @@ BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
     BOOST_CHECK_EQUAL( z[1], 0.8 );
     BOOST_CHECK_EQUAL( z[2], 0.4 );
 
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTL1 ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTL1 ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTL1 ), 0.2 + 0.2 + 0.4 );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTL1 ), 0.2 + 0.2 + 0.4 );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTLINF ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTLINF ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTLINF ), 0.4 );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTLINF ), 0.4 );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTTV ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTTV ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTKL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTKL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), INFINITY );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), INFINITY );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTHEL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTHEL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTL1 ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTL1 ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTL1 ), 0.2 + 0.2 + 0.4 );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTL1 ), 0.2 + 0.2 + 0.4 );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTLINF ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTLINF ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTLINF ), 0.4 );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTLINF ), 0.4 );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTTV ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTTV ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTKL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTKL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTKL ), INFINITY );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTKL ), INFINITY );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTHEL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTHEL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
     x.set( 1, 0.7 ); x.set( 2, 0.1 );
     y.set( 0, 0.1 ); y.set( 1, 0.5 );
-    BOOST_CHECK_CLOSE( dist( x, y, Prob::DISTKL ), 0.2 * std::log(0.2 / 0.1) + 0.7 * std::log(0.7 / 0.5) + 0.1 * std::log(0.1 / 0.4), tol );
-    BOOST_CHECK_CLOSE( dist( y, x, Prob::DISTKL ), 0.1 * std::log(0.1 / 0.2) + 0.5 * std::log(0.5 / 0.7) + 0.4 * std::log(0.4 / 0.1), tol );
+    BOOST_CHECK_CLOSE( dist( x, y, DISTKL ), 0.2 * std::log(0.2 / 0.1) + 0.7 * std::log(0.7 / 0.5) + 0.1 * std::log(0.1 / 0.4), tol );
+    BOOST_CHECK_CLOSE( dist( y, x, DISTKL ), 0.1 * std::log(0.1 / 0.2) + 0.5 * std::log(0.5 / 0.7) + 0.4 * std::log(0.4 / 0.1), tol );
 
     std::stringstream ss;
     ss << x;
index 8dcf61a..5df21ab 100644 (file)
@@ -32,18 +32,18 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     // check constructors
     Prob x1;
     BOOST_CHECK_EQUAL( x1.size(), 0 );
-    BOOST_CHECK( x1.p() == std::vector<Real>() );
+    BOOST_CHECK( x1.p() == Prob::container_type() );
 
     Prob x2( 3 );
     BOOST_CHECK_EQUAL( x2.size(), 3 );
-    BOOST_CHECK( x2.p() == std::vector<Real>( 3, 1.0 / 3.0 ) );
+    BOOST_CHECK( x2.p() == Prob::container_type( 3, 1.0 / 3.0 ) );
     BOOST_CHECK_EQUAL( x2[0], 1.0 / 3.0 );
     BOOST_CHECK_EQUAL( x2[1], 1.0 / 3.0 );
     BOOST_CHECK_EQUAL( x2[2], 1.0 / 3.0 );
 
     Prob x3( 4, 1.0 );
     BOOST_CHECK_EQUAL( x3.size(), 4 );
-    BOOST_CHECK( x3.p() == std::vector<Real>( 4, 1.0 ) );
+    BOOST_CHECK( x3.p() == Prob::container_type( 4, 1.0 ) );
     BOOST_CHECK_EQUAL( x3[0], 1.0 );
     BOOST_CHECK_EQUAL( x3[1], 1.0 );
     BOOST_CHECK_EQUAL( x3[2], 1.0 );
@@ -53,27 +53,32 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
     x3.set( 2, 2.0 );
     x3.set( 3, 4.0 );
 
-    Prob x4( x3.begin(), x3.end() );
+    std::vector<Real> v;
+    v.push_back( 0.5 );
+    v.push_back( 1.0 );
+    v.push_back( 2.0 );
+    v.push_back( 4.0 );
+    Prob x4( v.begin(), v.end(), 0 );
     BOOST_CHECK_EQUAL( x4.size(), 4 );
     BOOST_CHECK( x4.p() == x3.p() );
+    BOOST_CHECK( x4 == x3 );
     BOOST_CHECK_EQUAL( x4[0], 0.5 );
     BOOST_CHECK_EQUAL( x4[1], 1.0 );
     BOOST_CHECK_EQUAL( x4[2], 2.0 );
     BOOST_CHECK_EQUAL( x4[3], 4.0 );
 
-    x3.p() = std::vector<Real>( 4, 2.5 );
-    Prob x5( x3.begin(), x3.end(), x3.size() );
+    Prob x5( v.begin(), v.end(), v.size() );
     BOOST_CHECK_EQUAL( x5.size(), 4 );
     BOOST_CHECK( x5.p() == x3.p() );
-    BOOST_CHECK_EQUAL( x5[0], 2.5 );
-    BOOST_CHECK_EQUAL( x5[1], 2.5 );
-    BOOST_CHECK_EQUAL( x5[2], 2.5 );
-    BOOST_CHECK_EQUAL( x5[3], 2.5 );
+    BOOST_CHECK( x5 == x3 );
+    BOOST_CHECK_EQUAL( x5[0], 0.5 );
+    BOOST_CHECK_EQUAL( x5[1], 1.0 );
+    BOOST_CHECK_EQUAL( x5[2], 2.0 );
+    BOOST_CHECK_EQUAL( x5[3], 4.0 );
 
     std::vector<int> y( 3, 2 );
     Prob x6( y );
     BOOST_CHECK_EQUAL( x6.size(), 3 );
-    BOOST_CHECK( x6.p() == std::vector<Real>( 3, 2.0 ) );
     BOOST_CHECK_EQUAL( x6[0], 2.0 );
     BOOST_CHECK_EQUAL( x6[1], 2.0 );
     BOOST_CHECK_EQUAL( x6[2], 2.0 );
@@ -86,6 +91,7 @@ BOOST_AUTO_TEST_CASE( ConstructorsTest ) {
 }
 
 
+#ifndef DAI_SPARSE
 BOOST_AUTO_TEST_CASE( IteratorTest ) {
     Prob x( 5, 0.0 );
     size_t i;
@@ -116,6 +122,7 @@ BOOST_AUTO_TEST_CASE( IteratorTest ) {
     for( Prob::const_reverse_iterator crit = x.rbegin(); crit != x.rend(); crit++, i++ )
         BOOST_CHECK_EQUAL( *crit, 2 * i );
 }
+#endif
 
 
 BOOST_AUTO_TEST_CASE( QueriesTest ) {
@@ -125,37 +132,37 @@ BOOST_AUTO_TEST_CASE( QueriesTest ) {
 
     // test accumulate, min, max, sum, sumAbs, maxAbs
     BOOST_CHECK_EQUAL( x.sum(), 0.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_id<Real>() ), 0.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 1.0, std::plus<Real>(), fo_id<Real>() ), 1.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -1.0, std::plus<Real>(), fo_id<Real>() ), -1.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( 0.0, fo_id<Real>() ), 0.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( 1.0, fo_id<Real>() ), 1.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( -1.0, fo_id<Real>() ), -1.0 );
     BOOST_CHECK_EQUAL( x.max(), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -INFINITY, fo_max<Real>(), fo_id<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_id<Real>() ), 3.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -5.0, fo_max<Real>(), fo_id<Real>() ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -INFINITY, fo_id<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 3.0, fo_id<Real>(), false ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -5.0, fo_id<Real>(), false ), 2.0 );
     BOOST_CHECK_EQUAL( x.min(), -2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( INFINITY, fo_min<Real>(), fo_id<Real>() ), -2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_min<Real>(), fo_id<Real>() ), -3.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 5.0, fo_min<Real>(), fo_id<Real>() ), -2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( INFINITY, fo_id<Real>(), true ), -2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -3.0, fo_id<Real>(), true ), -3.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 5.0, fo_id<Real>(), true ), -2.0 );
     BOOST_CHECK_EQUAL( x.sumAbs(), 6.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_abs<Real>() ), 6.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 1.0, std::plus<Real>(), fo_abs<Real>() ), 7.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -1.0, std::plus<Real>(), fo_abs<Real>() ), 7.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( 0.0, fo_abs<Real>() ), 6.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( 1.0, fo_abs<Real>() ), 7.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( -1.0, fo_abs<Real>() ), 7.0 );
     BOOST_CHECK_EQUAL( x.maxAbs(), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 0.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 0.0, fo_abs<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 1.0, fo_abs<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -1.0, fo_abs<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 3.0, fo_abs<Real>(), false ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -3.0, fo_abs<Real>(), false ), 3.0 );
     x.set( 1, 1.0 );
     BOOST_CHECK_EQUAL( x.maxAbs(), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 0.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -1.0, fo_max<Real>(), fo_abs<Real>() ), 2.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( -3.0, fo_max<Real>(), fo_abs<Real>() ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 0.0, fo_abs<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 1.0, fo_abs<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -1.0, fo_abs<Real>(), false ), 2.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( 3.0, fo_abs<Real>(), false ), 3.0 );
+    BOOST_CHECK_EQUAL( x.accumulateMax( -3.0, fo_abs<Real>(), false ), 3.0 );
     for( size_t i = 0; i < x.size(); i++ )
         x.set( i, i ? (1.0 / i) : 0.0 );
-    BOOST_CHECK_EQUAL( x.accumulate( 0.0, std::plus<Real>(), fo_inv0<Real>() ), 10.0 );
+    BOOST_CHECK_EQUAL( x.accumulateSum( 0.0, fo_inv0<Real>() ), 10.0 );
     x /= x.sum();
 
     // test entropy
@@ -316,13 +323,13 @@ BOOST_AUTO_TEST_CASE( UnaryTransformationsTest ) {
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 0.5 );
 
-    y = x.normalized( Prob::NORMPROB );
+    y = x.normalized( NORMPROB );
     BOOST_CHECK_EQUAL( y[0], 0.5 );
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 0.5 );
 
     x.set( 0, -2.0 );
-    y = x.normalized( Prob::NORMLINF );
+    y = x.normalized( NORMLINF );
     BOOST_CHECK_EQUAL( y[0], -1.0 );
     BOOST_CHECK_EQUAL( y[1], 0.0 );
     BOOST_CHECK_EQUAL( y[2], 1.0 );
@@ -379,14 +386,14 @@ BOOST_AUTO_TEST_CASE( UnaryOperationsTest ) {
     BOOST_CHECK( x == y );
 
     x = xorg;
-    BOOST_CHECK_EQUAL( x.normalize( Prob::NORMPROB ), 3.0 );
+    BOOST_CHECK_EQUAL( x.normalize( NORMPROB ), 3.0 );
     BOOST_CHECK( x == y );
 
     y.set( 0, 2.0 / 2.0 );
     y.set( 1, 0.0 / 2.0 );
     y.set( 2, 1.0 / 2.0 );
     x = xorg;
-    BOOST_CHECK_EQUAL( x.normalize( Prob::NORMLINF ), 2.0 );
+    BOOST_CHECK_EQUAL( x.normalize( NORMLINF ), 2.0 );
     BOOST_CHECK( x == y );
 
     xorg.set( 0, -2.0 );
@@ -686,52 +693,65 @@ BOOST_AUTO_TEST_CASE( RelatedFunctionsTest ) {
     BOOST_CHECK_EQUAL( z[1], 0.8 );
     BOOST_CHECK_EQUAL( z[2], 0.4 );
 
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTL1 ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTL1 ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTL1 ), 0.2 + 0.2 + 0.4 );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTL1 ), 0.2 + 0.2 + 0.4 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTL1 ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTL1 ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTLINF ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTLINF ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTLINF ), 0.4 );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTLINF ), 0.4 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTLINF ), x.innerProduct( y, 0.0, fo_max<Real>(), fo_absdiff<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTLINF ), y.innerProduct( x, 0.0, fo_max<Real>(), fo_absdiff<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTTV ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTTV ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTTV ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) / 2.0 );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTTV ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) / 2.0 );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTKL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTKL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), INFINITY );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), INFINITY );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( x, x, Prob::DISTHEL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( y, y, Prob::DISTHEL ), 0.0 );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTHEL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTHEL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTL1 ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTL1 ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTL1 ), 0.2 + 0.2 + 0.4 );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTL1 ), 0.2 + 0.2 + 0.4 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTL1 ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTL1 ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTLINF ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTLINF ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTLINF ), 0.4 );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTLINF ), 0.4 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTLINF ), x.innerProduct( y, 0.0, fo_max<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTLINF ), y.innerProduct( x, 0.0, fo_max<Real>(), fo_absdiff<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTTV ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTTV ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTTV ), 0.5 * (0.2 + 0.2 + 0.4) );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTTV ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTTV ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_absdiff<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTKL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTKL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTKL ), INFINITY );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTKL ), INFINITY );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTKL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTKL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( x, x, DISTHEL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( y, y, DISTHEL ), 0.0 );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTHEL ), 0.5 * (0.2 + std::pow(std::sqrt(0.8) - std::sqrt(0.6), 2.0) + 0.4) );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTHEL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTHEL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_Hellinger<Real>() ) / 2.0 );
     x.set( 1, 0.7 ); x.set( 2, 0.1 );
     y.set( 0, 0.1 ); y.set( 1, 0.5 );
-    BOOST_CHECK_CLOSE( dist( x, y, Prob::DISTKL ), 0.2 * std::log(0.2 / 0.1) + 0.7 * std::log(0.7 / 0.5) + 0.1 * std::log(0.1 / 0.4), tol );
-    BOOST_CHECK_CLOSE( dist( y, x, Prob::DISTKL ), 0.1 * std::log(0.1 / 0.2) + 0.5 * std::log(0.5 / 0.7) + 0.4 * std::log(0.4 / 0.1), tol );
-    BOOST_CHECK_EQUAL( dist( x, y, Prob::DISTKL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
-    BOOST_CHECK_EQUAL( dist( y, x, Prob::DISTKL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
-
+    BOOST_CHECK_CLOSE( dist( x, y, DISTKL ), 0.2 * std::log(0.2 / 0.1) + 0.7 * std::log(0.7 / 0.5) + 0.1 * std::log(0.1 / 0.4), tol );
+    BOOST_CHECK_CLOSE( dist( y, x, DISTKL ), 0.1 * std::log(0.1 / 0.2) + 0.5 * std::log(0.5 / 0.7) + 0.4 * std::log(0.4 / 0.1), tol );
+    BOOST_CHECK_EQUAL( dist( x, y, DISTKL ), x.innerProduct( y, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+    BOOST_CHECK_EQUAL( dist( y, x, DISTKL ), y.innerProduct( x, 0.0, std::plus<Real>(), fo_KL<Real>() ) );
+
+    Prob xx(4), yy(4);
+    for( size_t i = 0; i < 3; i++ ) {
+        xx.set( i, x[i] );
+        yy.set( i, y[i] );
+    }
     std::stringstream ss;
-    ss << x;
+    ss << xx;
     std::string s;
     std::getline( ss, s );
-    BOOST_CHECK_EQUAL( s, std::string("(0.2, 0.7, 0.1)") );
+#ifdef DAI_SPARSE
+    BOOST_CHECK_EQUAL( s, std::string("(size:4, def:0.25, 0:0.2, 1:0.7, 2:0.1)") );
+#else
+    BOOST_CHECK_EQUAL( s, std::string("(0.2, 0.7, 0.1, 0.25)") );
+#endif
     std::stringstream ss2;
-    ss2 << y;
+    ss2 << yy;
     std::getline( ss2, s );
-    BOOST_CHECK_EQUAL( s, std::string("(0.1, 0.5, 0.4)") );
+#ifdef DAI_SPARSE
+    BOOST_CHECK_EQUAL( s, std::string("(size:4, def:0.25, 0:0.1, 1:0.5, 2:0.4)") );
+#else
+    BOOST_CHECK_EQUAL( s, std::string("(0.1, 0.5, 0.4, 0.25)") );
+#endif
 
     z = min( x, y );
     BOOST_CHECK_EQUAL( z[0], 0.1 );