include/dai/prob.h

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  libDAI is licensed under the terms of the GNU General Public License version
   4  *  2, or (at your option) any later version. libDAI is distributed without any
   5  *  warranty. See the file COPYING for more details.
   6  *
   7  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
   8  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
   9  */
  10
  11
  12 /// \file
  13 /// \brief Defines TProb<> and Prob classes which represent (probability) vectors (e.g., probability distributions of discrete random variables)
  14
  15
  16 #ifndef __defined_libdai_prob_h
  17 #define __defined_libdai_prob_h
  18
  19
  20 #include <cmath>
  21 #include <vector>
  22 #include <ostream>
  23 #include <algorithm>
  24 #include <numeric>
  25 #include <functional>
  26 #include <dai/util.h>
  27 #include <dai/exceptions.h>
  28
  29
  30 namespace dai {
  31
  32
  33 /// Function object that returns the value itself
  34 template<typename T> struct fo_id : public std::unary_function<T, T> {
  35     /// Returns \a x
  36     T operator()( const T &x ) const {
  37         return x;
  38     }
  39 };
  40
  41
  42 /// Function object that takes the absolute value
  43 template<typename T> struct fo_abs : public std::unary_function<T, T> {
  44     /// Returns abs(\a x)
  45     T operator()( const T &x ) const {
  46         if( x < (T)0 )
  47             return -x;
  48         else
  49             return x;
  50     }
  51 };
  52
  53
  54 /// Function object that takes the exponent
  55 template<typename T> struct fo_exp : public std::unary_function<T, T> {
  56     /// Returns exp(\a x)
  57     T operator()( const T &x ) const {
  58         return exp( x );
  59     }
  60 };
  61
  62
  63 /// Function object that takes the logarithm
  64 template<typename T> struct fo_log : public std::unary_function<T, T> {
  65     /// Returns log(\a x)
  66     T operator()( const T &x ) const {
  67         return log( x );
  68     }
  69 };
  70
  71
  72 /// Function object that takes the logarithm, except that log(0) is defined to be 0
  73 template<typename T> struct fo_log0 : public std::unary_function<T, T> {
  74     /// Returns (\a x == 0 ? 0 : log(\a x))
  75     T operator()( const T &x ) const {
  76         if( x )
  77             return log( x );
  78         else
  79             return 0;
  80     }
  81 };
  82
  83
  84 /// Function object that takes the inverse
  85 template<typename T> struct fo_inv : public std::unary_function<T, T> {
  86     /// Returns 1 / \a x
  87     T operator()( const T &x ) const {
  88         return 1 / x;
  89     }
  90 };
  91
  92
  93 /// Function object that takes the inverse, except that 1/0 is defined to be 0
  94 template<typename T> struct fo_inv0 : public std::unary_function<T, T> {
  95     /// Returns (\a x == 0 ? 0 : (1 / \a x))
  96     T operator()( const T &x ) const {
  97         if( x )
  98             return 1 / x;
  99         else
 100             return 0;
 101     }
 102 };
 103
 104
 105 /// Function object that returns p*log0(p)
 106 template<typename T> struct fo_plog0p : public std::unary_function<T, T> {
 107     /// Returns \a p * log0(\a p)
 108     T operator()( const T &p ) const {
 109         return p * dai::log0(p);
 110     }
 111 };
 112
 113
 114 /// Function object similar to std::divides(), but different in that dividing by zero results in zero
 115 template<typename T> struct fo_divides0 : public std::binary_function<T, T, T> {
 116     /// Returns (\a y == 0 ? 0 : (\a x / \a y))
 117     T operator()( const T &x, const T &y ) const {
 118         if( y == (T)0 )
 119             return (T)0;
 120         else
 121             return x / y;
 122     }
 123 };
 124
 125
 126 /// Function object useful for calculating the KL distance
 127 template<typename T> struct fo_KL : public std::binary_function<T, T, T> {
 128     /// Returns (\a p == 0 ? 0 : (\a p * (log(\a p) - log(\a q))))
 129     T operator()( const T &p, const T &q ) const {
 130         if( p == (T)0 )
 131             return (T)0;
 132         else
 133             return p * (log(p) - log(q));
 134     }
 135 };
 136
 137
 138 /// Function object useful for calculating the Hellinger distance
 139 template<typename T> struct fo_Hellinger : public std::binary_function<T, T, T> {
 140     /// Returns (sqrt(\a p) - sqrt(\a q))^2
 141     T operator()( const T &p, const T &q ) const {
 142         T x = sqrt(p) - sqrt(q);
 143         return x * x;
 144     }
 145 };
 146
 147
 148 /// Function object that returns x to the power y
 149 template<typename T> struct fo_pow : public std::binary_function<T, T, T> {
 150     /// Returns (\a x ^ \a y)
 151     T operator()( const T &x, const T &y ) const {
 152         if( y != 1 )
 153             return std::pow( x, y );
 154         else
 155             return x;
 156     }
 157 };
 158
 159
 160 /// Function object that returns the maximum of two values
 161 template<typename T> struct fo_max : public std::binary_function<T, T, T> {
 162     /// Returns (\a x > y ? x : y)
 163     T operator()( const T &x, const T &y ) const {
 164         return (x > y) ? x : y;
 165     }
 166 };
 167
 168
 169 /// Function object that returns the minimum of two values
 170 template<typename T> struct fo_min : public std::binary_function<T, T, T> {
 171     /// Returns (\a x > y ? y : x)
 172     T operator()( const T &x, const T &y ) const {
 173         return (x > y) ? y : x;
 174     }
 175 };
 176
 177
 178 /// Function object that returns the absolute difference of x and y
 179 template<typename T> struct fo_absdiff : public std::binary_function<T, T, T> {
 180     /// Returns abs( \a x - \a y )
 181     T operator()( const T &x, const T &y ) const {
 182         return dai::abs( x - y );
 183     }
 184 };
 185
 186
 187 /// Represents a vector with entries of type \a T.
 188 /** It is simply a <tt>std::vector</tt><<em>T</em>> with an interface designed for dealing with probability mass functions.
 189  *
 190  *  It is mainly used for representing measures on a finite outcome space, for example, the probability
 191  *  distribution of a discrete random variable. However, entries are not necessarily non-negative; it is also used to
 192  *  represent logarithms of probability mass functions.
 193  *
 194  *  \tparam T Should be a scalar that is castable from and to dai::Real and should support elementary arithmetic operations.
 195  */
 196 template <typename T> class TProb {
 197     private:
 198         /// The vector
 199         std::vector<T> _p;
 200
 201     public:
 202         /// Enumerates different ways of normalizing a probability measure.
 203         /**
 204          *  - NORMPROB means that the sum of all entries should be 1;
 205          *  - NORMLINF means that the maximum absolute value of all entries should be 1.
 206          */
 207         typedef enum { NORMPROB, NORMLINF } NormType;
 208         /// Enumerates different distance measures between probability measures.
 209         /**
 210          *  - DISTL1 is the \f$\ell_1\f$ distance (sum of absolute values of pointwise difference);
 211          *  - DISTLINF is the \f$\ell_\infty\f$ distance (maximum absolute value of pointwise difference);
 212          *  - DISTTV is the total variation distance (half of the \f$\ell_1\f$ distance);
 213          *  - DISTKL is the Kullback-Leibler distance (\f$\sum_i p_i (\log p_i - \log q_i)\f$).
 214          *  - DISTHEL is the Hellinger distance (\f$\frac{1}{2}\sum_i (\sqrt{p_i}-\sqrt{q_i})^2\f$).
 215          */
 216         typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL, DISTHEL } DistType;
 217
 218     /// \name Constructors and destructors
 219     //@{
 220         /// Default constructor (constructs empty vector)
 221         TProb() : _p() {}
 222
 223         /// 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$)
 224         explicit TProb( size_t n ) : _p(std::vector<T>(n, (T)1 / n)) {}
 225
 226         /// Construct vector of length \a n with each entry set to \a p
 227         explicit TProb( size_t n, T p ) : _p(n, p) {}
 228
 229         /// Construct vector from a range
 230         /** \tparam TIterator Iterates over instances that can be cast to \a T
 231          *  \param begin Points to first instance to be added.
 232          *  \param end Points just beyond last instance to be added.
 233          *  \param sizeHint For efficiency, the number of entries can be speficied by \a sizeHint.
 234          */
 235         template <typename TIterator>
 236         TProb( TIterator begin, TIterator end, size_t sizeHint=0 ) : _p() {
 237             _p.reserve( sizeHint );
 238             _p.insert( _p.begin(), begin, end );
 239         }
 240
 241         /// Construct vector from another vector
 242         /** \tparam S type of elements in \a v (should be castable to type \a T)
 243          *  \param v vector used for initialization
 244          */
 245         template <typename S>
 246         TProb( const std::vector<S> &v ) : _p() {
 247             _p.reserve( v.size() );
 248             _p.insert( _p.begin(), v.begin(), v.end() );
 249         }
 250     //@}
 251
 252         /// Constant iterator over the elements
 253         typedef typename std::vector<T>::const_iterator const_iterator;
 254         /// Iterator over the elements
 255         typedef typename std::vector<T>::iterator iterator;
 256         /// Constant reverse iterator over the elements
 257         typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
 258         /// Reverse iterator over the elements
 259         typedef typename std::vector<T>::reverse_iterator reverse_iterator;
 260
 261     /// \name Iterator interface
 262     //@{
 263         /// Returns iterator that points to the first element
 264         iterator begin() { return _p.begin(); }
 265         /// Returns constant iterator that points to the first element
 266         const_iterator begin() const { return _p.begin(); }
 267
 268         /// Returns iterator that points beyond the last element
 269         iterator end() { return _p.end(); }
 270         /// Returns constant iterator that points beyond the last element
 271         const_iterator end() const { return _p.end(); }
 272
 273         /// Returns reverse iterator that points to the last element
 274         reverse_iterator rbegin() { return _p.rbegin(); }
 275         /// Returns constant reverse iterator that points to the last element
 276         const_reverse_iterator rbegin() const { return _p.rbegin(); }
 277
 278         /// Returns reverse iterator that points beyond the first element
 279         reverse_iterator rend() { return _p.rend(); }
 280         /// Returns constant reverse iterator that points beyond the first element
 281         const_reverse_iterator rend() const { return _p.rend(); }
 282     //@}
 283
 284     /// \name Queries
 285     //@{
 286         /// Returns a const reference to the wrapped vector
 287         const std::vector<T> & p() const { return _p; }
 288
 289         /// Returns a reference to the wrapped vector
 290         std::vector<T> & p() { return _p; }
 291
 292         /// Returns a copy of the \a i 'th entry
 293         T operator[]( size_t i ) const {
 294 #ifdef DAI_DEBUG
 295             return _p.at(i);
 296 #else
 297             return _p[i];
 298 #endif
 299         }
 300
 301         /// Returns reference to the \a i 'th entry
 302         T& operator[]( size_t i ) { return _p[i]; }
 303
 304         /// Returns length of the vector (i.e., the number of entries)
 305         size_t size() const { return _p.size(); }
 306
 307         /// Accumulate over all values, similar to std::accumulate
 308         template<typename binOp, typename unOp> T accumulate( T init, binOp op1, unOp op2 ) const {
 309             T t = init;
 310             for( const_iterator it = begin(); it != end(); it++ )
 311                 t = op1( t, op2(*it) );
 312             return t;
 313         }
 314
 315         /// Returns the Shannon entropy of \c *this, \f$-\sum_i p_i \log p_i\f$
 316         T entropy() const { return -accumulate( (T)0, std::plus<T>(), fo_plog0p<T>() ); }
 317
 318         /// Returns maximum value of all entries
 319         T max() const { return accumulate( (T)(-INFINITY), fo_max<T>(), fo_id<T>() ); }
 320
 321         /// Returns minimum value of all entries
 322         T min() const { return accumulate( (T)INFINITY, fo_min<T>(), fo_id<T>() ); }
 323
 324         /// Returns sum of all entries
 325         T sum() const { return accumulate( (T)0, std::plus<T>(), fo_id<T>() ); }
 326
 327         /// Return sum of absolute value of all entries
 328         T sumAbs() const { return accumulate( (T)0, std::plus<T>(), fo_abs<T>() ); }
 329
 330         /// Returns maximum absolute value of all entries
 331         T maxAbs() const { return accumulate( (T)0, fo_max<T>(), fo_abs<T>() ); }
 332
 333         /// Returns \c true if one or more entries are NaN
 334         bool hasNaNs() const {
 335             bool foundnan = false;
 336             for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
 337                 if( isnan( *x ) ) {
 338                     foundnan = true;
 339                     break;
 340                 }
 341             return foundnan;
 342         }
 343
 344         /// Returns \c true if one or more entries are negative
 345         bool hasNegatives() const {
 346             return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<T>(), (T)0 ) ) != _p.end());
 347         }
 348
 349         /// Returns a pair consisting of the index of the maximum value and the maximum value itself
 350         std::pair<size_t,T> argmax() const {
 351             T max = _p[0];
 352             size_t arg = 0;
 353             for( size_t i = 1; i < size(); i++ ) {
 354               if( _p[i] > max ) {
 355                 max = _p[i];
 356                 arg = i;
 357               }
 358             }
 359             return std::make_pair(arg,max);
 360         }
 361
 362         /// Returns a random index, according to the (normalized) distribution described by *this
 363         size_t draw() {
 364             Real x = rnd_uniform() * sum();
 365             T s = 0;
 366             for( size_t i = 0; i < size(); i++ ) {
 367                 s += _p[i];
 368                 if( s > x )
 369                     return i;
 370             }
 371             return( size() - 1 );
 372         }
 373
 374         /// Lexicographical comparison
 375         /** \pre <tt>this->size() == q.size()</tt>
 376          */
 377         bool operator<= (const TProb<T> & q) const {
 378             DAI_DEBASSERT( size() == q.size() );
 379             return lexicographical_compare( begin(), end(), q.begin(), q.end() );
 380         }
 381     //@}
 382
 383     /// \name Unary transformations
 384     //@{
 385         /// Returns the result of applying operation \a op pointwise on \c *this
 386         template<typename unaryOp> TProb<T> pwUnaryTr( unaryOp op ) const {
 387             TProb<T> r;
 388             r._p.reserve( size() );
 389             std::transform( _p.begin(), _p.end(), back_inserter( r._p ), op );
 390             return r;
 391         }
 392
 393         /// Returns pointwise absolute value
 394         TProb<T> abs() const { return pwUnaryTr( fo_abs<T>() ); }
 395
 396         /// Returns pointwise exponent
 397         TProb<T> exp() const { return pwUnaryTr( fo_exp<T>() ); }
 398
 399         /// Returns pointwise logarithm
 400         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
 401          */
 402         TProb<T> log(bool zero=false) const {
 403             if( zero )
 404                 return pwUnaryTr( fo_log0<T>() );
 405             else
 406                 return pwUnaryTr( fo_log<T>() );
 407         }
 408
 409         /// Returns pointwise inverse
 410         /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
 411          */
 412         TProb<T> inverse(bool zero=true) const {
 413             if( zero )
 414                 return pwUnaryTr( fo_inv0<T>() );
 415             else
 416                 return pwUnaryTr( fo_inv<T>() );
 417         }
 418
 419         /// Returns normalized copy of \c *this, using the specified norm
 420         /** \throw NOT_NORMALIZABLE if the norm is zero
 421          */
 422         TProb<T> normalized( NormType norm = NORMPROB ) const {
 423             T Z = 0;
 424             if( norm == NORMPROB )
 425                 Z = sum();
 426             else if( norm == NORMLINF )
 427                 Z = maxAbs();
 428             if( Z == (T)0 ) {
 429                 DAI_THROW(NOT_NORMALIZABLE);
 430                 return *this;
 431             } else
 432                 return pwUnaryTr( std::bind2nd( std::divides<T>(), Z ) );
 433         }
 434     //@}
 435
 436     /// \name Unary operations
 437     //@{
 438         /// Applies unary operation \a op pointwise
 439         template<typename unaryOp> TProb<T>& pwUnaryOp( unaryOp op ) {
 440             std::transform( _p.begin(), _p.end(), _p.begin(), op );
 441             return *this;
 442         }
 443
 444         /// Draws all entries i.i.d. from a uniform distribution on [0,1)
 445         TProb<T>& randomize() {
 446             std::generate( _p.begin(), _p.end(), rnd_uniform );
 447             return *this;
 448         }
 449
 450         /// Sets all entries to \f$1/n\f$ where \a n is the length of the vector
 451         TProb<T>& setUniform () {
 452             fill( (T)1 / size() );
 453             return *this;
 454         }
 455
 456         /// Applies absolute value pointwise
 457         const TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
 458
 459         /// Applies exponent pointwise
 460         const TProb<T>& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
 461
 462         /// Applies logarithm pointwise
 463         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
 464          */
 465         const TProb<T>& takeLog(bool zero=false) {
 466             if( zero ) {
 467                 return pwUnaryOp( fo_log0<T>() );
 468             } else
 469                 return pwUnaryOp( fo_log<T>() );
 470         }
 471
 472         /// Normalizes vector using the specified norm
 473         /** \throw NOT_NORMALIZABLE if the norm is zero
 474          */
 475         T normalize( NormType norm=NORMPROB ) {
 476             T Z = 0;
 477             if( norm == NORMPROB )
 478                 Z = sum();
 479             else if( norm == NORMLINF )
 480                 Z = maxAbs();
 481             if( Z == (T)0 )
 482                 DAI_THROW(NOT_NORMALIZABLE);
 483             else
 484                 *this /= Z;
 485             return Z;
 486         }
 487     //@}
 488
 489     /// \name Operations with scalars
 490     //@{
 491         /// Sets all entries to \a x
 492         TProb<T> & fill(T x) {
 493             std::fill( _p.begin(), _p.end(), x );
 494             return *this;
 495         }
 496
 497         /// Adds scalar \a x to each entry
 498         TProb<T>& operator+= (T x) {
 499             if( x != 0 )
 500                 return pwUnaryOp( std::bind2nd( std::plus<T>(), x ) );
 501             else
 502                 return *this;
 503         }
 504
 505         /// Subtracts scalar \a x from each entry
 506         TProb<T>& operator-= (T x) {
 507             if( x != 0 )
 508                 return pwUnaryOp( std::bind2nd( std::minus<T>(), x ) );
 509             else
 510                 return *this;
 511         }
 512
 513         /// Multiplies each entry with scalar \a x
 514         TProb<T>& operator*= (T x) {
 515             if( x != 1 )
 516                 return pwUnaryOp( std::bind2nd( std::multiplies<T>(), x ) );
 517             else
 518                 return *this;
 519         }
 520
 521         /// Divides each entry by scalar \a x
 522         TProb<T>& operator/= (T x) {
 523             DAI_DEBASSERT( x != 0 );
 524             if( x != 1 )
 525                 return pwUnaryOp( std::bind2nd( std::divides<T>(), x ) );
 526             else
 527                 return *this;
 528         }
 529
 530         /// Raises entries to the power \a x
 531         TProb<T>& operator^= (T x) {
 532             if( x != (T)1 )
 533                 return pwUnaryOp( std::bind2nd( fo_pow<T>(), x) );
 534             else
 535                 return *this;
 536         }
 537     //@}
 538
 539     /// \name Transformations with scalars
 540     //@{
 541         /// Returns sum of \c *this and scalar \a x
 542         TProb<T> operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }
 543
 544         /// Returns difference of \c *this and scalar \a x
 545         TProb<T> operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }
 546
 547         /// Returns product of \c *this with scalar \a x
 548         TProb<T> operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }
 549
 550         /// Returns quotient of \c *this and scalar \a x, where division by 0 yields 0
 551         TProb<T> operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }
 552
 553         /// Returns \c *this raised to the power \a x
 554         TProb<T> operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }
 555     //@}
 556
 557     /// \name Operations with other equally-sized vectors
 558     //@{
 559         /// Applies binary operation pointwise on two vectors
 560         /** \tparam binaryOp Type of function object that accepts two arguments of type \a T and outputs a type \a T
 561          *  \param q Right operand
 562          *  \param op Operation of type \a binaryOp
 563          */
 564         template<typename binaryOp> TProb<T>& pwBinaryOp( const TProb<T> &q, binaryOp op ) {
 565             DAI_DEBASSERT( size() == q.size() );
 566             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), op );
 567             return *this;
 568         }
 569
 570         /// Pointwise addition with \a q
 571         /** \pre <tt>this->size() == q.size()</tt>
 572          */
 573         TProb<T>& operator+= (const TProb<T> & q) { return pwBinaryOp( q, std::plus<T>() ); }
 574
 575         /// Pointwise subtraction of \a q
 576         /** \pre <tt>this->size() == q.size()</tt>
 577          */
 578         TProb<T>& operator-= (const TProb<T> & q) { return pwBinaryOp( q, std::minus<T>() ); }
 579
 580         /// Pointwise multiplication with \a q
 581         /** \pre <tt>this->size() == q.size()</tt>
 582          */
 583         TProb<T>& operator*= (const TProb<T> & q) { return pwBinaryOp( q, std::multiplies<T>() ); }
 584
 585         /// Pointwise division by \a q, where division by 0 yields 0
 586         /** \pre <tt>this->size() == q.size()</tt>
 587          *  \see divide(const TProb<T> &)
 588          */
 589         TProb<T>& operator/= (const TProb<T> & q) { return pwBinaryOp( q, fo_divides0<T>() ); }
 590
 591         /// Pointwise division by \a q, where division by 0 yields +Inf
 592         /** \pre <tt>this->size() == q.size()</tt>
 593          *  \see operator/=(const TProb<T> &)
 594          */
 595         TProb<T>& divide (const TProb<T> & q) { return pwBinaryOp( q, std::divides<T>() ); }
 596
 597         /// Pointwise power
 598         /** \pre <tt>this->size() == q.size()</tt>
 599          */
 600         TProb<T>& operator^= (const TProb<T> & q) { return pwBinaryOp( q, fo_pow<T>() ); }
 601     //@}
 602
 603     /// \name Transformations with other equally-sized vectors
 604     //@{
 605         /// Returns the result of applying binary operation \a op pointwise on \c *this and \a q
 606         /** \tparam binaryOp Type of function object that accepts two arguments of type \a T and outputs a type \a T
 607          *  \param q Right operand
 608          *  \param op Operation of type \a binaryOp
 609          */
 610         template<typename binaryOp> TProb<T> pwBinaryTr( const TProb<T> &q, binaryOp op ) const {
 611             DAI_DEBASSERT( size() == q.size() );
 612             TProb<T> r;
 613             r._p.reserve( size() );
 614             std::transform( _p.begin(), _p.end(), q._p.begin(), back_inserter( r._p ), op );
 615             return r;
 616         }
 617
 618         /// Returns sum of \c *this and \a q
 619         /** \pre <tt>this->size() == q.size()</tt>
 620          */
 621         TProb<T> operator+ ( const TProb<T>& q ) const { return pwBinaryTr( q, std::plus<T>() ); }
 622
 623         /// Return \c *this minus \a q
 624         /** \pre <tt>this->size() == q.size()</tt>
 625          */
 626         TProb<T> operator- ( const TProb<T>& q ) const { return pwBinaryTr( q, std::minus<T>() ); }
 627
 628         /// Return product of \c *this with \a q
 629         /** \pre <tt>this->size() == q.size()</tt>
 630          */
 631         TProb<T> operator* ( const TProb<T> &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }
 632
 633         /// Returns quotient of \c *this with \a q, where division by 0 yields 0
 634         /** \pre <tt>this->size() == q.size()</tt>
 635          *  \see divided_by(const TProb<T> &)
 636          */
 637         TProb<T> operator/ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }
 638
 639         /// Pointwise division by \a q, where division by 0 yields +Inf
 640         /** \pre <tt>this->size() == q.size()</tt>
 641          *  \see operator/(const TProb<T> &)
 642          */
 643         TProb<T> divided_by( const TProb<T> &q ) const { return pwBinaryTr( q, std::divides<T>() ); }
 644
 645         /// Returns \c *this to the power \a q
 646         /** \pre <tt>this->size() == q.size()</tt>
 647          */
 648         TProb<T> operator^ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }
 649     //@}
 650
 651         /// Performs a generalized inner product, similar to std::inner_product
 652         /** \pre <tt>this->size() == q.size()</tt>
 653          */
 654         template<typename binOp1, typename binOp2> T innerProduct( const TProb<T> &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
 655             DAI_DEBASSERT( size() == q.size() );
 656             return std::inner_product( begin(), end(), q.begin(), init, binaryOp1, binaryOp2 );
 657         }
 658 };
 659
 660
 661 /// Returns distance between \a p and \a q, measured using distance measure \a dt
 662 /** \relates TProb
 663  *  \pre <tt>this->size() == q.size()</tt>
 664  */
 665 template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
 666     switch( dt ) {
 667         case TProb<T>::DISTL1:
 668             return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() );
 669         case TProb<T>::DISTLINF:
 670             return p.innerProduct( q, (T)0, fo_max<T>(), fo_absdiff<T>() );
 671         case TProb<T>::DISTTV:
 672             return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() ) / 2;
 673         case TProb<T>::DISTKL:
 674             return p.innerProduct( q, (T)0, std::plus<T>(), fo_KL<T>() );
 675         case TProb<T>::DISTHEL:
 676             return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
 677         default:
 678             DAI_THROW(UNKNOWN_ENUM_VALUE);
 679             return INFINITY;
 680     }
 681 }
 682
 683
 684 /// Writes a TProb<T> to an output stream
 685 /** \relates TProb
 686  */
 687 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
 688     os << "[";
 689     std::copy( p.p().begin(), p.p().end(), std::ostream_iterator<T>(os, " ") );
 690     os << "]";
 691     return os;
 692 }
 693
 694
 695 /// Returns the pointwise minimum of \a a and \a b
 696 /** \relates TProb
 697  *  \pre <tt>this->size() == q.size()</tt>
 698  */
 699 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
 700     return a.pwBinaryTr( b, fo_min<T>() );
 701 }
 702
 703
 704 /// Returns the pointwise maximum of \a a and \a b
 705 /** \relates TProb
 706  *  \pre <tt>this->size() == q.size()</tt>
 707  */
 708 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
 709     return a.pwBinaryTr( b, fo_max<T>() );
 710 }
 711
 712
 713 /// Represents a vector with entries of type dai::Real.
 714 typedef TProb<Real> Prob;
 715
 716
 717 } // end of namespace dai
 718
 719
 720 #endif