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 negative of \c *this
 394         TProb<T> operator- () const { return pwUnaryTr( std::negate<T>() ); }
 395
 396         /// Returns pointwise absolute value
 397         TProb<T> abs() const { return pwUnaryTr( fo_abs<T>() ); }
 398
 399         /// Returns pointwise exponent
 400         TProb<T> exp() const { return pwUnaryTr( fo_exp<T>() ); }
 401
 402         /// Returns pointwise logarithm
 403         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
 404          */
 405         TProb<T> log(bool zero=false) const {
 406             if( zero )
 407                 return pwUnaryTr( fo_log0<T>() );
 408             else
 409                 return pwUnaryTr( fo_log<T>() );
 410         }
 411
 412         /// Returns pointwise inverse
 413         /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
 414          */
 415         TProb<T> inverse(bool zero=true) const {
 416             if( zero )
 417                 return pwUnaryTr( fo_inv0<T>() );
 418             else
 419                 return pwUnaryTr( fo_inv<T>() );
 420         }
 421
 422         /// Returns normalized copy of \c *this, using the specified norm
 423         /** \throw NOT_NORMALIZABLE if the norm is zero
 424          */
 425         TProb<T> normalized( NormType norm = NORMPROB ) const {
 426             T Z = 0;
 427             if( norm == NORMPROB )
 428                 Z = sum();
 429             else if( norm == NORMLINF )
 430                 Z = maxAbs();
 431             if( Z == (T)0 ) {
 432                 DAI_THROW(NOT_NORMALIZABLE);
 433                 return *this;
 434             } else
 435                 return pwUnaryTr( std::bind2nd( std::divides<T>(), Z ) );
 436         }
 437     //@}
 438
 439     /// \name Unary operations
 440     //@{
 441         /// Applies unary operation \a op pointwise
 442         template<typename unaryOp> TProb<T>& pwUnaryOp( unaryOp op ) {
 443             std::transform( _p.begin(), _p.end(), _p.begin(), op );
 444             return *this;
 445         }
 446
 447         /// Draws all entries i.i.d. from a uniform distribution on [0,1)
 448         TProb<T>& randomize() {
 449             std::generate( _p.begin(), _p.end(), rnd_uniform );
 450             return *this;
 451         }
 452
 453         /// Sets all entries to \f$1/n\f$ where \a n is the length of the vector
 454         TProb<T>& setUniform () {
 455             fill( (T)1 / size() );
 456             return *this;
 457         }
 458
 459         /// Applies absolute value pointwise
 460         const TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
 461
 462         /// Applies exponent pointwise
 463         const TProb<T>& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
 464
 465         /// Applies logarithm pointwise
 466         /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
 467          */
 468         const TProb<T>& takeLog(bool zero=false) {
 469             if( zero ) {
 470                 return pwUnaryOp( fo_log0<T>() );
 471             } else
 472                 return pwUnaryOp( fo_log<T>() );
 473         }
 474
 475         /// Normalizes vector using the specified norm
 476         /** \throw NOT_NORMALIZABLE if the norm is zero
 477          */
 478         T normalize( NormType norm=NORMPROB ) {
 479             T Z = 0;
 480             if( norm == NORMPROB )
 481                 Z = sum();
 482             else if( norm == NORMLINF )
 483                 Z = maxAbs();
 484             if( Z == (T)0 )
 485                 DAI_THROW(NOT_NORMALIZABLE);
 486             else
 487                 *this /= Z;
 488             return Z;
 489         }
 490     //@}
 491
 492     /// \name Operations with scalars
 493     //@{
 494         /// Sets all entries to \a x
 495         TProb<T> & fill(T x) {
 496             std::fill( _p.begin(), _p.end(), x );
 497             return *this;
 498         }
 499
 500         /// Adds scalar \a x to each entry
 501         TProb<T>& operator+= (T x) {
 502             if( x != 0 )
 503                 return pwUnaryOp( std::bind2nd( std::plus<T>(), x ) );
 504             else
 505                 return *this;
 506         }
 507
 508         /// Subtracts scalar \a x from each entry
 509         TProb<T>& operator-= (T x) {
 510             if( x != 0 )
 511                 return pwUnaryOp( std::bind2nd( std::minus<T>(), x ) );
 512             else
 513                 return *this;
 514         }
 515
 516         /// Multiplies each entry with scalar \a x
 517         TProb<T>& operator*= (T x) {
 518             if( x != 1 )
 519                 return pwUnaryOp( std::bind2nd( std::multiplies<T>(), x ) );
 520             else
 521                 return *this;
 522         }
 523
 524         /// Divides each entry by scalar \a x
 525         TProb<T>& operator/= (T x) {
 526             DAI_DEBASSERT( x != 0 );
 527             if( x != 1 )
 528                 return pwUnaryOp( std::bind2nd( std::divides<T>(), x ) );
 529             else
 530                 return *this;
 531         }
 532
 533         /// Raises entries to the power \a x
 534         TProb<T>& operator^= (T x) {
 535             if( x != (T)1 )
 536                 return pwUnaryOp( std::bind2nd( fo_pow<T>(), x) );
 537             else
 538                 return *this;
 539         }
 540     //@}
 541
 542     /// \name Transformations with scalars
 543     //@{
 544         /// Returns sum of \c *this and scalar \a x
 545         TProb<T> operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }
 546
 547         /// Returns difference of \c *this and scalar \a x
 548         TProb<T> operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }
 549
 550         /// Returns product of \c *this with scalar \a x
 551         TProb<T> operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }
 552
 553         /// Returns quotient of \c *this and scalar \a x, where division by 0 yields 0
 554         TProb<T> operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }
 555
 556         /// Returns \c *this raised to the power \a x
 557         TProb<T> operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }
 558     //@}
 559
 560     /// \name Operations with other equally-sized vectors
 561     //@{
 562         /// Applies binary operation pointwise on two vectors
 563         /** \tparam binaryOp Type of function object that accepts two arguments of type \a T and outputs a type \a T
 564          *  \param q Right operand
 565          *  \param op Operation of type \a binaryOp
 566          */
 567         template<typename binaryOp> TProb<T>& pwBinaryOp( const TProb<T> &q, binaryOp op ) {
 568             DAI_DEBASSERT( size() == q.size() );
 569             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), op );
 570             return *this;
 571         }
 572
 573         /// Pointwise addition with \a q
 574         /** \pre <tt>this->size() == q.size()</tt>
 575          */
 576         TProb<T>& operator+= (const TProb<T> & q) { return pwBinaryOp( q, std::plus<T>() ); }
 577
 578         /// Pointwise subtraction of \a q
 579         /** \pre <tt>this->size() == q.size()</tt>
 580          */
 581         TProb<T>& operator-= (const TProb<T> & q) { return pwBinaryOp( q, std::minus<T>() ); }
 582
 583         /// Pointwise multiplication with \a q
 584         /** \pre <tt>this->size() == q.size()</tt>
 585          */
 586         TProb<T>& operator*= (const TProb<T> & q) { return pwBinaryOp( q, std::multiplies<T>() ); }
 587
 588         /// Pointwise division by \a q, where division by 0 yields 0
 589         /** \pre <tt>this->size() == q.size()</tt>
 590          *  \see divide(const TProb<T> &)
 591          */
 592         TProb<T>& operator/= (const TProb<T> & q) { return pwBinaryOp( q, fo_divides0<T>() ); }
 593
 594         /// Pointwise division by \a q, where division by 0 yields +Inf
 595         /** \pre <tt>this->size() == q.size()</tt>
 596          *  \see operator/=(const TProb<T> &)
 597          */
 598         TProb<T>& divide (const TProb<T> & q) { return pwBinaryOp( q, std::divides<T>() ); }
 599
 600         /// Pointwise power
 601         /** \pre <tt>this->size() == q.size()</tt>
 602          */
 603         TProb<T>& operator^= (const TProb<T> & q) { return pwBinaryOp( q, fo_pow<T>() ); }
 604     //@}
 605
 606     /// \name Transformations with other equally-sized vectors
 607     //@{
 608         /// Returns the result of applying binary operation \a op pointwise on \c *this and \a q
 609         /** \tparam binaryOp Type of function object that accepts two arguments of type \a T and outputs a type \a T
 610          *  \param q Right operand
 611          *  \param op Operation of type \a binaryOp
 612          */
 613         template<typename binaryOp> TProb<T> pwBinaryTr( const TProb<T> &q, binaryOp op ) const {
 614             DAI_DEBASSERT( size() == q.size() );
 615             TProb<T> r;
 616             r._p.reserve( size() );
 617             std::transform( _p.begin(), _p.end(), q._p.begin(), back_inserter( r._p ), op );
 618             return r;
 619         }
 620
 621         /// Returns sum of \c *this and \a q
 622         /** \pre <tt>this->size() == q.size()</tt>
 623          */
 624         TProb<T> operator+ ( const TProb<T>& q ) const { return pwBinaryTr( q, std::plus<T>() ); }
 625
 626         /// Return \c *this minus \a q
 627         /** \pre <tt>this->size() == q.size()</tt>
 628          */
 629         TProb<T> operator- ( const TProb<T>& q ) const { return pwBinaryTr( q, std::minus<T>() ); }
 630
 631         /// Return product of \c *this with \a q
 632         /** \pre <tt>this->size() == q.size()</tt>
 633          */
 634         TProb<T> operator* ( const TProb<T> &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }
 635
 636         /// Returns quotient of \c *this with \a q, where division by 0 yields 0
 637         /** \pre <tt>this->size() == q.size()</tt>
 638          *  \see divided_by(const TProb<T> &)
 639          */
 640         TProb<T> operator/ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }
 641
 642         /// Pointwise division by \a q, where division by 0 yields +Inf
 643         /** \pre <tt>this->size() == q.size()</tt>
 644          *  \see operator/(const TProb<T> &)
 645          */
 646         TProb<T> divided_by( const TProb<T> &q ) const { return pwBinaryTr( q, std::divides<T>() ); }
 647
 648         /// Returns \c *this to the power \a q
 649         /** \pre <tt>this->size() == q.size()</tt>
 650          */
 651         TProb<T> operator^ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }
 652     //@}
 653
 654         /// Performs a generalized inner product, similar to std::inner_product
 655         /** \pre <tt>this->size() == q.size()</tt>
 656          */
 657         template<typename binOp1, typename binOp2> T innerProduct( const TProb<T> &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
 658             DAI_DEBASSERT( size() == q.size() );
 659             return std::inner_product( begin(), end(), q.begin(), init, binaryOp1, binaryOp2 );
 660         }
 661 };
 662
 663
 664 /// Returns distance between \a p and \a q, measured using distance measure \a dt
 665 /** \relates TProb
 666  *  \pre <tt>this->size() == q.size()</tt>
 667  */
 668 template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
 669     switch( dt ) {
 670         case TProb<T>::DISTL1:
 671             return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() );
 672         case TProb<T>::DISTLINF:
 673             return p.innerProduct( q, (T)0, fo_max<T>(), fo_absdiff<T>() );
 674         case TProb<T>::DISTTV:
 675             return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() ) / 2;
 676         case TProb<T>::DISTKL:
 677             return p.innerProduct( q, (T)0, std::plus<T>(), fo_KL<T>() );
 678         case TProb<T>::DISTHEL:
 679             return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
 680         default:
 681             DAI_THROW(UNKNOWN_ENUM_VALUE);
 682             return INFINITY;
 683     }
 684 }
 685
 686
 687 /// Writes a TProb<T> to an output stream
 688 /** \relates TProb
 689  */
 690 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
 691     os << "[";
 692     std::copy( p.p().begin(), p.p().end(), std::ostream_iterator<T>(os, " ") );
 693     os << "]";
 694     return os;
 695 }
 696
 697
 698 /// Returns the pointwise minimum of \a a and \a b
 699 /** \relates TProb
 700  *  \pre <tt>this->size() == q.size()</tt>
 701  */
 702 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
 703     return a.pwBinaryTr( b, fo_min<T>() );
 704 }
 705
 706
 707 /// Returns the pointwise maximum of \a a and \a b
 708 /** \relates TProb
 709  *  \pre <tt>this->size() == q.size()</tt>
 710  */
 711 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
 712     return a.pwBinaryTr( b, fo_max<T>() );
 713 }
 714
 715
 716 /// Represents a vector with entries of type dai::Real.
 717 typedef TProb<Real> Prob;
 718
 719
 720 } // end of namespace dai
 721
 722
 723 #endif