Fixed testem failure caused by rounding error

[libdai.git] / include / dai / prob.h

diff --git a/include/dai/prob.h b/include/dai/prob.h
index 1291242..4b85374 100644 (file)
--- a/include/dai/prob.h
+++ b/include/dai/prob.h
@@ -1,6 +1,7 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
-    Radboud University Nijmegen, The Netherlands
-    
+/*  Copyright (C) 2006-2008  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
+    Radboud University Nijmegen, The Netherlands /
+    Max Planck Institute for Biological Cybernetics, Germany
+

     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify


@@ -19,160 +20,229 @@
 */
 
 
 */
 
 

+/// \file
+/// \brief Defines TProb<T> and Prob classes
+/// \todo Rename to Vector<T>
+
+

 #ifndef __defined_libdai_prob_h
 #define __defined_libdai_prob_h
 
 
 #ifndef __defined_libdai_prob_h
 #define __defined_libdai_prob_h
 
 

-#include <complex>

 #include <cmath>
 #include <vector>
 #include <cmath>
 #include <vector>

-#include <iostream>
+#include <ostream>

 #include <cassert>
 #include <cassert>

+#include <algorithm>
+#include <numeric>
+#include <functional>

 #include <dai/util.h>
 #include <dai/util.h>

+#include <dai/exceptions.h>

 
 
 namespace dai {
 
 
 
 
 namespace dai {
 
 

-typedef double                  Real;
-typedef std::complex<double>    Complex;
-
-template<typename T> class      TProb;
-typedef TProb<Real>             Prob;
-typedef TProb<Complex>          CProb;
-
-
-/// TProb<T> implements a probability vector of type T.
-/// T should be castable from and to double and to complex.
+/// 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.
+ */

 template <typename T> class TProb {
 template <typename T> class TProb {

-    protected:
-        /// The entries
-        std::vector<T> _p;
-

     private:
     private:

-        /// Calculate x times log(x), or 0 if x == 0
-        Complex xlogx( Real x ) const { return( x == 0.0 ? 0.0 : Complex(x) * std::log(Complex(x))); }
+        /// The vector
+        std::vector<T> _p;

 
     public:
 
     public:

-        /// NORMPROB means that the sum of all entries should be 1
-        /// NORMLINF means that the maximum absolute value of all entries should be 1
+        /// 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;
+         *  - NORMLINF means that the maximum absolute value of all entries should be 1.
+         */

         typedef enum { NORMPROB, NORMLINF } NormType;
         typedef enum { NORMPROB, NORMLINF } NormType;

-        /// 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
-        typedef enum { DISTL1, DISTLINF, DISTTV } DistType;
+        /// 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.
+         */
+        typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;

         
         /// Default constructor
         TProb() : _p() {}
         
         
         /// Default constructor
         TProb() : _p() {}
         

-        /// Construct uniform distribution of given length
-        TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
+        /// 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 with given length and initial value
-        TProb( size_t n, Real p ) : _p(std::vector<T>(n,(T)p)) {}
+        /// Construct vector of length n with each entry set to p
+        explicit TProb( size_t n, Real p ) : _p(n, (T)p) {}

         
         

-        /// Construct with given length and initial array
-        TProb( size_t n, const Real* p ) {
-            // Reserve-push_back is faster than resize-copy
-            _p.reserve( n );
-            for( size_t i = 0; i < n; i++ ) 
-                _p.push_back( p[i] );
+        /// Construct vector from a range
+        /** \tparam Iterator Iterates over instances that can be cast to T.
+         *  \param begin Points to first instance to be added.
+         *  \param end Points just beyond last instance to be added.
+         *  \param sizeHint For efficiency, the number of entries can be speficied by sizeHint.
+         */
+        template <typename Iterator>
+        TProb( Iterator begin, Iterator end, size_t sizeHint=0 ) : _p() {
+            _p.reserve( sizeHint );
+            _p.insert( _p.begin(), begin, end );

         }
         
         }
         

-        /// Copy constructor
-        TProb( const TProb<T> & x ) : _p(x._p) {}
-        
-        /// Assignment operator
-        TProb<T> & operator=( const TProb<T> &x ) {
-            if( this != &x ) {
-                _p = x._p;
-            }
-            return *this;
-        }
-        
-        /// Provide read access to _p
+        /// Returns a const reference to the vector

         const std::vector<T> & p() const { return _p; }
 
         const std::vector<T> & p() const { return _p; }
 

-        /// Provide full access to _p
+        /// Returns a reference to the vector

         std::vector<T> & p() { return _p; }
         
         std::vector<T> & p() { return _p; }
         

-        /// Provide read access to ith element of _p
-        T operator[]( size_t i ) const { return _p[i]; }
+        /// Returns a copy of the i'th entry
+        T operator[]( size_t i ) const { 
+#ifdef DAI_DEBUG
+            return _p.at(i);
+#else
+            return _p[i];
+#endif
+        }

         
         

-        /// Provide full access to ith element of _p
+        /// Returns reference to the i'th entry

         T& operator[]( size_t i ) { return _p[i]; }
         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(); }

 
 

-        /// Set all elements to x
+        /// Returns iterator pointing beyond last entry
+        iterator end() { return _p.end(); }
+
+        /// Returns const iterator pointing beyond last entry
+        const_iterator end() const { return _p.end(); }
+
+        /// Sets all entries to x

         TProb<T> & fill(T x) { 
         TProb<T> & fill(T x) { 

-            for( size_t i = 0; i < size(); i++ )
-                _p[i] = x;
+            std::fill( _p.begin(), _p.end(), x );

             return *this;
         }
 
             return *this;
         }
 

-        /// Set all elements to iid random numbers from uniform(0,1) distribution
+        /// Draws all entries i.i.d. from a uniform distribution on [0,1)

         TProb<T> & randomize() { 
         TProb<T> & randomize() { 

-            for( size_t i = 0; i < size(); i++ )
-                _p[i] = rnd_uniform();
+            std::generate(_p.begin(), _p.end(), rnd_uniform);

             return *this;
         }
 
             return *this;
         }
 

-        /// Return size
+        /// Returns length of the vector, i.e., the number of entries

         size_t size() const {
             return _p.size();
         }
 
         size_t size() const {
             return _p.size();
         }
 

-        /// Make entries zero if (Real) absolute value smaller than epsilon
-        TProb<T>& makeZero (Real epsilon) {
+        /// Sets entries that are smaller than epsilon to 0
+        TProb<T>& makeZero( Real epsilon ) {

             for( size_t i = 0; i < size(); i++ )
             for( size_t i = 0; i < size(); i++ )

-                if( fabs((Real)_p[i]) < epsilon )
+                if( fabs(_p[i]) < epsilon )

                     _p[i] = 0;
             return *this;
         }
                     _p[i] = 0;
             return *this;
         }

+        
+        /// Set all entries to 1.0/size()
+        TProb<T>& setUniform () {
+            fill(1.0/size());
+            return *this;
+        }

 
 

-        /// Multiplication with T x
-        TProb<T>& operator*= (T x) {
+        /// Sets entries that are smaller than epsilon to epsilon
+        TProb<T>& makePositive( Real epsilon ) {

             for( size_t i = 0; i < size(); i++ )
             for( size_t i = 0; i < size(); i++ )

-                _p[i] *= x;
+                if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
+                    _p[i] = epsilon;
+            return *this;
+        }
+
+        /// 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;
         }
 
             return *this;
         }
 

-        /// Return product of *this with T x
+        /// Returns product of *this with scalar x

         TProb<T> operator* (T x) const {
             TProb<T> prod( *this );
             prod *= x;
             return prod;
         }
 
         TProb<T> operator* (T x) const {
             TProb<T> prod( *this );
             prod *= x;
             return prod;
         }
 

-        /// Division by T x
+        /// Divides each entry by scalar x

         TProb<T>& operator/= (T x) {
         TProb<T>& operator/= (T x) {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( x != 0.0 );
 #endif
             assert( x != 0.0 );
 #endif

-            for( size_t i = 0; i < size(); i++ )
-                _p[i] /= x;
+            std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );

             return *this;
         }
 
             return *this;
         }
 

-        /// Return quotient of *this and T x
+        /// Returns quotient of *this and scalar x

         TProb<T> operator/ (T x) const {
         TProb<T> operator/ (T x) const {

-            TProb<T> prod( *this );
-            prod /= x;
-            return prod;
+            TProb<T> quot( *this );
+            quot /= x;
+            return quot;

         }
 
         }
 

-        /// Pointwise multiplication with q
-        TProb<T>& operator*= (const TProb<T> & q) {
-#ifdef DEBUG
+        /// 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 sum of *this and scalar x
+        TProb<T> operator+ (T x) const {
+            TProb<T> sum( *this );
+            sum += x;
+            return sum;
+        }
+
+        /// 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 difference of *this and scalar x
+        TProb<T> operator- (T x) const {
+            TProb<T> diff( *this );
+            diff -= x;
+            return diff;
+        }
+
+        /// Lexicographical comparison (sizes should be identical)
+        bool operator<= (const TProb<T> & q) const {
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             for( size_t i = 0; i < size(); i++ )
             assert( size() == q.size() );
 #endif
             for( size_t i = 0; i < size(); i++ )

-                _p[i] *= q[i];
+                if( !(_p[i] <= q[i]) )
+                    return false;
+            return true;
+        }
+
+        /// Pointwise multiplication with q (sizes should be identical)
+        TProb<T>& operator*= (const TProb<T> & q) {
+#ifdef DAI_DEBUG
+            assert( size() == q.size() );
+#endif
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );

             return *this;
         }
         
             return *this;
         }
         

-        /// Return product of *this with q
+        /// Return product of *this with q (sizes should be identical)

         TProb<T> operator* (const TProb<T> & q) const {
         TProb<T> operator* (const TProb<T> & q) const {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             TProb<T> prod( *this );
             assert( size() == q.size() );
 #endif
             TProb<T> prod( *this );


@@ -180,29 +250,18 @@ template <typename T> class TProb {
             return prod;
         }
 
             return prod;
         }
 

-        /// Pointwise addition with q
+        /// Pointwise addition with q (sizes should be identical)

         TProb<T>& operator+= (const TProb<T> & q) {
         TProb<T>& operator+= (const TProb<T> & q) {

-#ifdef DEBUG
-            assert( size() == q.size() );
-#endif
-            for( size_t i = 0; i < size(); i++ )
-                _p[i] += q[i];
-            return *this;
-        }
-        
-        /// Pointwise subtraction of q
-        TProb<T>& operator-= (const TProb<T> & q) {
-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             assert( size() == q.size() );
 #endif

-            for( size_t i = 0; i < size(); i++ )
-                _p[i] -= q[i];
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );

             return *this;
         }
         
             return *this;
         }
         

-        /// Return sum of *this and q
+        /// Returns sum of *this and q (sizes should be identical)

         TProb<T> operator+ (const TProb<T> & q) const {
         TProb<T> operator+ (const TProb<T> & q) const {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             TProb<T> sum( *this );
             assert( size() == q.size() );
 #endif
             TProb<T> sum( *this );


@@ -210,26 +269,32 @@ template <typename T> class TProb {
             return sum;
         }
         
             return sum;
         }
         

-        /// Return *this minus q
+        /// Pointwise subtraction of q (sizes should be identical)
+        TProb<T>& operator-= (const TProb<T> & q) {
+#ifdef DAI_DEBUG
+            assert( size() == q.size() );
+#endif
+            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 {
         TProb<T> operator- (const TProb<T> & q) const {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             assert( size() == q.size() );
 #endif

-            TProb<T> sum( *this );
-            sum -= q;
-            return sum;
+            TProb<T> diff( *this );
+            diff -= q;
+            return diff;

         }
 
         }
 

-        /// Pointwise division by q
+        /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)

         TProb<T>& operator/= (const TProb<T> & q) {
         TProb<T>& operator/= (const TProb<T> & q) {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             for( size_t i = 0; i < size(); i++ ) {
             assert( size() == q.size() );
 #endif
             for( size_t i = 0; i < size(); i++ ) {

-#ifdef DEBUG
-//              assert( q[i] != 0.0 );
-#endif
-                if( q[i] == 0.0 )       // FIXME
+                if( q[i] == 0.0 )

                     _p[i] = 0.0;
                 else
                     _p[i] /= q[i];
                     _p[i] = 0.0;
                 else
                     _p[i] /= q[i];


@@ -237,19 +302,18 @@ template <typename T> class TProb {
             return *this;
         }
         
             return *this;
         }
         

-        /// Pointwise division by q, division by zero yields infinity
+        /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)

         TProb<T>& divide (const TProb<T> & q) {
         TProb<T>& divide (const TProb<T> & q) {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             assert( size() == q.size() );
 #endif

-            for( size_t i = 0; i < size(); i++ )
-                _p[i] /= q[i];
+            std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );

             return *this;
         }
         
             return *this;
         }
         

-        /// Return quotient of *this with q
+        /// Returns quotient of *this with q (sizes should be identical)

         TProb<T> operator/ (const TProb<T> & q) const {
         TProb<T> operator/ (const TProb<T> & q) const {

-#ifdef DEBUG
+#ifdef DAI_DEBUG

             assert( size() == q.size() );
 #endif
             TProb<T> quot( *this );
             assert( size() == q.size() );
 #endif
             TProb<T> quot( *this );


@@ -257,136 +321,110 @@ template <typename T> class TProb {
             return quot;
         }
 
             return quot;
         }
 

-        /// Return pointwise inverse
-        TProb<T> inverse(bool zero = false) const {
+        /// 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
             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++ ) {
-#ifdef DEBUG
-                    assert( _p[i] != 0.0 );
-#endif
+                for( size_t i = 0; i < size(); i++ )

                     inv._p.push_back( 1.0 / _p[i] );
                     inv._p.push_back( 1.0 / _p[i] );

-                }

             return inv;
         }
 
             return inv;
         }
 

-        /// Return *this to the power of a (pointwise)
+        /// Raises entries to the power a

         TProb<T>& operator^= (Real a) {
         TProb<T>& operator^= (Real a) {

-            if( a != 1.0 ) {
-                for( size_t i = 0; i < size(); i++ )
-                    _p[i] = std::pow( _p[i], 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;
         }
 
             return *this;
         }
 

-        /// Pointwise power of a
+        /// Returns *this raised to the power a

         TProb<T> operator^ (Real a) const {
         TProb<T> operator^ (Real a) const {

-            TProb<T> power;
-            if( a != 1.0 ) {
-                power._p.reserve( size() );
-                for( size_t i = 0; i < size(); i++ )
-                    power._p.push_back( std::pow( _p[i], a ) );
-            } else
-                power = *this;
+            TProb<T> power(*this);
+            power ^= a;

             return power;
         }
 
             return power;
         }
 

-        /// Pointwise exp
-        TProb<T> exp() const {
-            TProb<T> e;
-            e._p.reserve( size() );
-            for( size_t i = 0; i < size(); i++ )
-                e._p.push_back( std::exp( _p[i] ) );
-            return e;
+        /// Returns pointwise signum
+        TProb<T> sgn() const {
+            TProb<T> x;
+            x._p.reserve( size() );
+            for( size_t i = 0; i < size(); i++ ) {
+                T s = 0;
+                if( _p[i] > 0 )
+                    s = 1;
+                else if( _p[i] < 0 )
+                    s = -1;
+                x._p.push_back( s );
+            }
+            return x;

         }
 
         }
 

-        /// Pointwise log
-        TProb<T> log() const {
-            TProb<T> l;
-            l._p.reserve( size() );
+        /// Returns pointwise absolute value
+        TProb<T> abs() const {
+            TProb<T> x;
+            x._p.reserve( size() );

             for( size_t i = 0; i < size(); i++ )
             for( size_t i = 0; i < size(); i++ )

-                l._p.push_back( std::log( _p[i] ) );
-            return l;
+                x._p.push_back( _p[i] < 0 ? (-_p[i]) : _p[i] );
+            return x;

         }
 
         }
 

-        /// Pointwise log (or 0 if == 0)
-        TProb<T> log0() const {
-            TProb<T> l0;
-            l0._p.reserve( size() );
-            for( size_t i = 0; i < size(); i++ )
-                l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( _p[i] ) );
-            return l0;
+        /// Applies exp pointwise
+        const TProb<T>& takeExp() {
+            std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::exp) );
+            return *this;

         }
 
         }
 

-        /// Pointwise (complex) log (or 0 if == 0)
-/*      CProb clog0() const {
-            CProb l0;
-            l0._p.reserve( size() );
-            for( size_t i = 0; i < size(); i++ )
-                l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( Complex( _p[i] ) ) );
-            return l0;
-        }*/
-
-        /// Return distance of p and q
-        friend Real dist( const TProb<T> & p, const TProb<T> & q, DistType dt ) {
-#ifdef DEBUG
-            assert( p.size() == q.size() );
-#endif
-            Real result = 0.0;
-            switch( dt ) {
-                case DISTL1:
-                    for( size_t i = 0; i < p.size(); i++ )
-                        result += fabs((Real)p[i] - (Real)q[i]);
-                    break;
-                    
-                case DISTLINF:
-                    for( size_t i = 0; i < p.size(); i++ ) {
-                        Real z = fabs((Real)p[i] - (Real)q[i]);
-                        if( z > result )
-                            result = z;
-                    }
-                    break;
+        /// 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;
+        }

 
 

-                case DISTTV:
-                    for( size_t i = 0; i < p.size(); i++ )
-                        result += fabs((Real)p[i] - (Real)q[i]);
-                    result *= 0.5;
-                    break;
-            }
-            return result;
+        /// Returns pointwise exp
+        TProb<T> exp() const {
+            TProb<T> e(*this);
+            e.takeExp();
+            return e;

         }
 
         }
 

-        /// Return (complex) Kullback-Leibler distance with q
-        friend Complex KL_dist( const TProb<T> & p, const TProb<T> & q ) {
-#ifdef DEBUG
-            assert( p.size() == q.size() );
-#endif
-            Complex result = 0.0;
-            for( size_t i = 0; i < p.size(); i++ ) {
-                if( (Real) p[i] != 0.0 ) {
-                    Complex p_i = p[i];
-                    Complex q_i = q[i];
-                    result += p_i * (std::log(p_i) - std::log(q_i));
-                }
-            }
-            return result;
+        /// Returns pointwise log
+        /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
+         */
+        TProb<T> log(bool zero=false) const {
+            TProb<T> l(*this);
+            l.takeLog(zero);
+            return l;

         }
 
         }
 

-        /// Return sum of all entries
-        T totalSum() const {
-            T Z = 0.0;
-            for( size_t i = 0; i < size(); i++ )
-                Z += _p[i];
+        /// Returns sum of all entries
+        T sum() const {
+            T Z = std::accumulate( _p.begin(),  _p.end(), (T)0 );

             return Z;
         }
 
             return Z;
         }
 

-        /// Converts entries to Real and returns maximum absolute value
+        /// Return sum of absolute value of all entries
+        T sumAbs() const {
+            T s = 0;
+            for( size_t i = 0; i < size(); i++ )
+                s += fabs( (Real) _p[i] );
+            return s;
+        }
+
+        /// Returns maximum absolute value of all entries

         T maxAbs() const {
         T maxAbs() const {

-            T Z = 0.0;
+            T Z = 0;

             for( size_t i = 0; i < size(); i++ ) {
                 Real mag = fabs( (Real) _p[i] );
                 if( mag > Z )
             for( size_t i = 0; i < size(); i++ ) {
                 Real mag = fabs( (Real) _p[i] );
                 if( mag > Z )


@@ -395,86 +433,173 @@ template <typename T> class TProb {
             return Z;
         }
 
             return Z;
         }
 

-        /// Returns maximum value
+        /// Returns maximum value of all entries

         T max() const {
         T max() const {

-            T Z = 0.0;
-            for( size_t i = 0; i < size(); i++ ) {
-                if( _p[i] > Z )
-                    Z = _p[i];
-            }
+            T Z = *std::max_element( _p.begin(), _p.end() );

             return Z;
         }
 
             return Z;
         }
 

-        /// Normalize, using the specified norm
-        T normalize( NormType norm ) {
+        /// Returns minimum value of all entries
+        T min() const {
+            T Z = *std::min_element( _p.begin(), _p.end() );
+            return Z;
+        }
+
+        /// Returns {arg,}maximum value
+        std::pair<size_t,T> argmax() const {
+            T max = _p[0];
+            size_t arg = 0;
+            for( size_t i = 1; i < size(); i++ ) {
+              if( _p[i] > max ) {
+                max = _p[i];
+                arg = i;
+              }
+            }
+            return std::make_pair(arg,max);
+        }
+
+        /// Normalizes vector using the specified norm
+        T normalize( NormType norm=NORMPROB ) {

             T Z = 0.0;
             if( norm == NORMPROB )
             T Z = 0.0;
             if( norm == NORMPROB )

-                Z = totalSum();
+                Z = sum();

             else if( norm == NORMLINF )
                 Z = maxAbs();
             else if( norm == NORMLINF )
                 Z = maxAbs();

-#ifdef DEBUG
-            assert( Z != 0.0 );
-#endif
-            T Zi = 1.0 / Z;
-            for( size_t i = 0; i < size(); i++ )
-               _p[i] *= Zi;
+            if( Z == 0.0 )
+                DAI_THROW(NOT_NORMALIZABLE);
+            else
+                *this /= Z;

             return Z;
         }
 
             return Z;
         }
 

-        /// Return normalized copy of *this, using the specified norm
-        TProb<T> normalized( NormType norm ) const {
-            T Z = 0.0;
-            if( norm == NORMPROB ) 
-                Z = totalSum();
-            else if( norm == NORMLINF )
-                Z = maxAbs();
-#ifdef DEBUG
-            assert( Z != 0.0 );
-#endif
-            Z = 1.0 / Z;
-            
-            TProb<T> result;
-            result._p.reserve( size() );
-            for( size_t i = 0; i < size(); i++ )
-                result._p.push_back( _p[i] * 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;
         }
     
         /// Returns true if one or more entries are NaN
         bool hasNaNs() const {
             return result;
         }
     
         /// Returns true if one or more entries are NaN
         bool hasNaNs() const {

-            bool NaNs = false;
-            for( size_t i = 0; i < size() && !NaNs; i++ ) 
-                if( isnan( _p[i] ) )
-                    NaNs = true;
-            return NaNs;
+            bool foundnan = false;
+            for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
+                if( isnan( *x ) ) {
+                    foundnan = true;
+                    break;
+                }
+            return foundnan;

         }
 
         /// Returns true if one or more entries are negative
         bool hasNegatives() const {
         }
 
         /// Returns true if one or more entries are negative
         bool hasNegatives() const {

-            bool Negatives = false;
-            for( size_t i = 0; i < size() && !Negatives; i++ ) 
-                if( _p[i] < 0.0 )
-                    Negatives = true;
-            return Negatives;
+            return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());

         }
         }

-
-        /// Returns (complex) entropy
-        Complex entropy() const {
-            Complex S = 0.0;
+        
+        /// Returns entropy of *this
+        Real entropy() const {
+            Real S = 0.0;

             for( size_t i = 0; i < size(); i++ )
             for( size_t i = 0; i < size(); i++ )

-                S -= xlogx(_p[i]);
+                S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));

             return S;
         }
 
             return S;
         }
 

-        friend std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
-            for( size_t i = 0; i < P.size(); i++ )
-                os << P._p[i] << " ";
-            os << std::endl;
-            return os;
+        /// Returns a random index, according to the (normalized) distribution described by *this
+        size_t draw() {
+            double 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 );

         }
 };
 
 
         }
 };
 
 

+/// Returns distance of p and q (sizes should be identical), measured using distance measure dt
+/** \relates TProb
+ */
+template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
+#ifdef DAI_DEBUG
+    assert( p.size() == q.size() );
+#endif
+    Real result = 0.0;
+    switch( dt ) {
+        case TProb<T>::DISTL1:
+            for( size_t i = 0; i < p.size(); i++ )
+                result += fabs((Real)p[i] - (Real)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]);
+                if( z > result )
+                    result = z;
+            }
+            break;
+
+        case TProb<T>::DISTTV:
+            for( size_t i = 0; i < p.size(); i++ )
+                result += fabs((Real)p[i] - (Real)q[i]);
+            result *= 0.5;
+            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]));
+            }
+    }
+    return result;
+}
+
+
+/// Writes a TProb<T> to an output stream
+/** \relates TProb
+ */
+template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
+    os << "[";
+    std::copy( P.p().begin(), P.p().end(), std::ostream_iterator<T>(os, " ") );
+    os << "]";
+    return os;
+}
+
+
+/// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
+/** \relates TProb
+ */
+template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
+    assert( a.size() == b.size() );
+    TProb<T> result( a.size() );
+    for( size_t i = 0; i < a.size(); i++ )
+        if( a[i] < b[i] )
+            result[i] = a[i];
+        else
+            result[i] = b[i];
+    return result;
+}
+
+
+/// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
+/** \relates TProb
+ */
+template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
+    assert( a.size() == b.size() );
+    TProb<T> result( a.size() );
+    for( size_t i = 0; i < a.size(); i++ )
+        if( a[i] > b[i] )
+            result[i] = a[i];
+        else
+            result[i] = b[i];
+    return result;
+}
+
+
+/// Represents a vector with entries of type Real.
+typedef TProb<Real> Prob;
+
+

 } // end of namespace dai
 
 
 } // end of namespace dai