Various cleanups and documentation improvements for SmallSet and Prob
[libdai.git] / include / dai / prob.h
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;