Extended SWIG python interface (inspired by Kyle Ellrott): inference is possible...
[libdai.git] / include / dai / prob.h
index d382478..51dfe47 100644 (file)
@@ -1,16 +1,13 @@
 /*  This file is part of libDAI - http://www.libdai.org/
  *
- *  libDAI is licensed under the terms of the GNU General Public License version
- *  2, or (at your option) any later version. libDAI is distributed without any
- *  warranty. See the file COPYING for more details.
+ *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
  *
- *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
- *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
  */
 
 
 /// \file
-/// \brief Defines TProb<> and Prob classes which represent (probability) vectors
+/// \brief Defines TProb<> and Prob classes which represent (probability) vectors (e.g., probability distributions of discrete random variables)
 
 
 #ifndef __defined_libdai_prob_h
@@ -135,12 +132,22 @@ template<typename T> struct fo_KL : public std::binary_function<T, T, T> {
 };
 
 
+/// Function object useful for calculating the Hellinger distance
+template<typename T> struct fo_Hellinger : public std::binary_function<T, T, T> {
+    /// Returns (sqrt(\a p) - sqrt(\a q))^2
+    T operator()( const T &p, const T &q ) const {
+        T x = sqrt(p) - sqrt(q);
+        return x * x;
+    }
+};
+
+
 /// Function object that returns x to the power y
 template<typename T> struct fo_pow : public std::binary_function<T, T, T> {
     /// Returns (\a x ^ \a y)
     T operator()( const T &x, const T &y ) const {
         if( y != 1 )
-            return std::pow( x, y );
+            return pow( x, y );
         else
             return x;
     }
@@ -183,53 +190,47 @@ template<typename T> struct fo_absdiff : public std::binary_function<T, T, T> {
  *
  *  \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 {
+template <typename T>
+class TProb {
+    public:
+        /// Type of data structure used for storing the values
+        typedef std::vector<T> container_type;
+
+        /// Shorthand
+        typedef TProb<T> this_type;
+
     private:
-        /// The vector
-        std::vector<T> _p;
+        /// The data structure that stores the values
+        container_type _p;
 
     public:
-        /// 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;
-        /// 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$).
-         */
-        typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
-
     /// \name Constructors and destructors
     //@{
         /// Default constructor (constructs empty vector)
         TProb() : _p() {}
 
         /// 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)) {}
+        explicit TProb( size_t n ) : _p( n, (T)1 / n ) {}
 
         /// Construct vector of length \a n with each entry set to \a p
-        explicit TProb( size_t n, T p ) : _p(n, p) {}
+        explicit TProb( size_t n, T p ) : _p( n, p ) {}
 
         /// Construct vector from a range
         /** \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.
          */
         template <typename TIterator>
-        TProb( TIterator begin, TIterator end, size_t sizeHint=0 ) : _p() {
+        TProb( TIterator begin, TIterator end, size_t sizeHint ) : _p() {
             _p.reserve( sizeHint );
             _p.insert( _p.begin(), begin, end );
         }
 
         /// 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
+         *  \param v vector used for initialization.
          */
         template <typename S>
         TProb( const std::vector<S> &v ) : _p() {
@@ -239,13 +240,13 @@ template <typename T> class TProb {
     //@}
 
         /// Constant iterator over the elements
-        typedef typename std::vector<T>::const_iterator const_iterator;
+        typedef typename container_type::const_iterator const_iterator;
         /// Iterator over the elements
-        typedef typename std::vector<T>::iterator iterator;
+        typedef typename container_type::iterator iterator;
         /// Constant reverse iterator over the elements
-        typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
+        typedef typename container_type::const_reverse_iterator const_reverse_iterator;
         /// Reverse iterator over the elements
-        typedef typename std::vector<T>::reverse_iterator reverse_iterator;
+        typedef typename container_type::reverse_iterator reverse_iterator;
 
     /// \name Iterator interface
     //@{
@@ -270,16 +271,17 @@ template <typename T> class TProb {
         const_reverse_iterator rend() const { return _p.rend(); }
     //@}
 
-    /// \name Queries
+    /// \name Miscellaneous operations
     //@{
-        /// Returns a const reference to the wrapped vector
-        const std::vector<T> & p() const { return _p; }
-
-        /// Returns a reference to the wrapped vector
-        std::vector<T> & p() { return _p; }
+        void resize( size_t sz ) {
+            _p.resize( sz );
+        }
+    //@}
 
-        /// Returns a copy of the \a i 'th entry
-        T operator[]( size_t i ) const {
+    /// \name Get/set individual entries
+    //@{
+        /// Gets \a i 'th entry
+        T get( size_t i ) const { 
 #ifdef DAI_DEBUG
             return _p.at(i);
 #else
@@ -287,43 +289,87 @@ template <typename T> class TProb {
 #endif
         }
 
-        /// Returns reference to the \a i 'th entry
-        T& operator[]( size_t i ) { return _p[i]; }
+        /// Sets \a i 'th entry to \a val
+        void set( size_t i, T val ) {
+            DAI_DEBASSERT( i < _p.size() );
+            _p[i] = val;
+        }
+    //@}
+
+    /// \name Queries
+    //@{
+        /// Returns a const reference to the wrapped container
+        const container_type& p() const { return _p; }
+
+        /// Returns a reference to the wrapped container
+        container_type& p() { return _p; }
+
+        /// Returns a copy of the \a i 'th entry
+        T operator[]( size_t i ) const { return get(i); }
 
         /// Returns length of the vector (i.e., the number of entries)
         size_t size() const { return _p.size(); }
 
-        /// Accumulate over all values, similar to std::accumulate
-        template<typename binOp, typename unOp> T accumulate( T init, binOp op1, unOp op2 ) const {
-            T t = init;
+        /// 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 = op1( t, op2(*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 {
             bool foundnan = false;
-            for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
-                if( isnan( *x ) ) {
+            for( const_iterator x = _p.begin(); x != _p.end(); x++ )
+                if( dai::isnan( *x ) ) {
                     foundnan = true;
                     break;
                 }
@@ -345,7 +391,7 @@ template <typename T> class TProb {
                 arg = i;
               }
             }
-            return std::make_pair(arg,max);
+            return std::make_pair( arg, max );
         }
 
         /// Returns a random index, according to the (normalized) distribution described by *this
@@ -353,7 +399,7 @@ template <typename T> class TProb {
             Real x = rnd_uniform() * sum();
             T s = 0;
             for( size_t i = 0; i < size(); i++ ) {
-                s += _p[i];
+                s += get(i);
                 if( s > x )
                     return i;
             }
@@ -363,50 +409,49 @@ template <typename T> 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 this_type& q ) const {
+            if( size() != q.size() )
+                return false;
+            return p() == q.p();
+        }
+
+        /// Formats a TProb as a string
+        std::string toString() const {
+            std::stringstream ss;
+            ss << *this;
+            return ss.str();
+        }
     //@}
 
     /// \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;
         }
 
-        // OBSOLETE
-        /// Returns pointwise signum
-        /** \deprecated Functionality was not widely used
-         */
-        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;
-        }
+        /// Returns negative of \c *this
+        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
@@ -416,7 +461,7 @@ template <typename T> 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
@@ -424,11 +469,13 @@ template <typename T> class TProb {
         }
 
         /// Returns normalized copy of \c *this, using the specified norm
-        TProb<T> normalized( NormType norm = NORMPROB ) const {
+        /** \throw NOT_NORMALIZABLE if the norm is zero
+         */
+        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);
@@ -441,33 +488,33 @@ template <typename T> 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
-        const TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
+        this_type& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
 
         /// Applies exponent pointwise
-        const 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>.
          */
-        const TProb<T>& takeLog(bool zero=false) {
+        this_type& takeLog(bool zero=false) {
             if( zero ) {
                 return pwUnaryOp( fo_log0<T>() );
             } else
@@ -475,11 +522,13 @@ template <typename T> class TProb {
         }
 
         /// Normalizes vector using the specified norm
-        T normalize( NormType norm=NORMPROB ) {
+        /** \throw NOT_NORMALIZABLE if the norm is zero
+         */
+        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);
@@ -492,35 +541,13 @@ template <typename T> 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;
         }
 
-        // OBSOLETE
-        /// Sets entries that are smaller (in absolute value) than \a epsilon to 0
-        /** \deprecated Functionality was not widely used
-         */
-        TProb<T>& makeZero( T epsilon ) {
-            for( size_t i = 0; i < size(); i++ )
-                if( (_p[i] < epsilon) && (_p[i] > -epsilon) )
-                    _p[i] = 0;
-            return *this;
-        }
-        
-        // OBSOLETE
-        /// Sets entries that are smaller than \a epsilon to \a epsilon
-        /** \deprecated Functionality was not widely used
-         */
-        TProb<T>& makePositive( T epsilon ) {
-            for( size_t i = 0; i < size(); i++ )
-                if( (0 < _p[i]) && (_p[i] < epsilon) )
-                    _p[i] = epsilon;
-            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
@@ -528,7 +555,7 @@ template <typename T> 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
@@ -536,24 +563,23 @@ template <typename T> 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
                 return *this;
         }
 
-        /// Divides each entry by scalar \a x
-        TProb<T>& operator/= (T x) {
-            DAI_DEBASSERT( x != 0 );
+        /// Divides each entry by scalar \a x, where division by 0 yields 0
+        this_type& operator/= (T x) {
             if( x != 1 )
-                return pwUnaryOp( std::bind2nd( std::divides<T>(), x ) );
+                return pwUnaryOp( std::bind2nd( fo_divides0<T>(), x ) );
             else
                 return *this;
         }
 
         /// 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
@@ -564,19 +590,19 @@ template <typename T> 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
@@ -586,7 +612,7 @@ template <typename T> 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;
@@ -595,34 +621,34 @@ template <typename T> 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
@@ -632,7 +658,7 @@ template <typename T> 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() );
@@ -643,40 +669,40 @@ template <typename T> 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 );
         }
@@ -687,16 +713,18 @@ template <typename T> 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 DISTHEL:
+            return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
         default:
             DAI_THROW(UNKNOWN_ENUM_VALUE);
             return INFINITY;
@@ -708,9 +736,10 @@ template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TPro
 /** \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 << "]";
+    os << "(";
+    for( size_t i = 0; i < p.size(); i++ )
+        os << ((i != 0) ? ", " : "") << p.get(i);
+    os << ")";
     return os;
 }