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<T> and Prob classes
  14 /// \todo Rename to Vector<T>
  15
  16
  17 #ifndef __defined_libdai_prob_h
  18 #define __defined_libdai_prob_h
  19
  20
  21 #include <cmath>
  22 #include <vector>
  23 #include <ostream>
  24 #include <cassert>
  25 #include <algorithm>
  26 #include <numeric>
  27 #include <functional>
  28 #include <dai/util.h>
  29 #include <dai/exceptions.h>
  30
  31
  32 namespace dai {
  33
  34
  35 /// Represents a vector with entries of type \a T.
  36 /** A TProb<T> is a std::vector<T> with an interface designed for dealing with probability mass functions.
  37  *  It is mainly used for representing measures on a finite outcome space, e.g., the probability
  38  *  distribution of a discrete random variable.
  39  *  \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
  40  */
  41 template <typename T> class TProb {
  42     private:
  43         /// The vector
  44         std::vector<T> _p;
  45
  46     public:
  47         /// Iterator over entries
  48         typedef typename std::vector<T>::iterator iterator;
  49         /// Const iterator over entries
  50         typedef typename std::vector<T>::const_iterator const_iterator;
  51
  52         /// Enumerates different ways of normalizing a probability measure.
  53         /**
  54          *  - NORMPROB means that the sum of all entries should be 1;
  55          *  - NORMLINF means that the maximum absolute value of all entries should be 1.
  56          */
  57         typedef enum { NORMPROB, NORMLINF } NormType;
  58         /// Enumerates different distance measures between probability measures.
  59         /**
  60          *  - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
  61          *  - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
  62          *  - DISTTV is the Total Variation distance;
  63          *  - DISTKL is the Kullback-Leibler distance.
  64          */
  65         typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
  66
  67         /// Default constructor
  68         TProb() : _p() {}
  69
  70         /// Construct uniform distribution over n outcomes, i.e., a vector of length n with each entry set to 1/n
  71         explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
  72
  73         /// Construct vector of length n with each entry set to p
  74         explicit TProb( size_t n, Real p ) : _p(n, (T)p) {}
  75
  76         /// Construct vector from a range
  77         /** \tparam Iterator Iterates over instances that can be cast to T.
  78          *  \param begin Points to first instance to be added.
  79          *  \param end Points just beyond last instance to be added.
  80          *  \param sizeHint For efficiency, the number of entries can be speficied by sizeHint.
  81          */
  82         template <typename Iterator>
  83         TProb( Iterator begin, Iterator end, size_t sizeHint=0 ) : _p() {
  84             _p.reserve( sizeHint );
  85             _p.insert( _p.begin(), begin, end );
  86         }
  87
  88         /// Returns a const reference to the vector
  89         const std::vector<T> & p() const { return _p; }
  90
  91         /// Returns a reference to the vector
  92         std::vector<T> & p() { return _p; }
  93
  94         /// Returns a copy of the i'th entry
  95         T operator[]( size_t i ) const {
  96 #ifdef DAI_DEBUG
  97             return _p.at(i);
  98 #else
  99             return _p[i];
 100 #endif
 101         }
 102
 103         /// Returns reference to the i'th entry
 104         T& operator[]( size_t i ) { return _p[i]; }
 105
 106         /// Returns iterator pointing to first entry
 107         iterator begin() { return _p.begin(); }
 108
 109         /// Returns const iterator pointing to first entry
 110         const_iterator begin() const { return _p.begin(); }
 111
 112         /// Returns iterator pointing beyond last entry
 113         iterator end() { return _p.end(); }
 114
 115         /// Returns const iterator pointing beyond last entry
 116         const_iterator end() const { return _p.end(); }
 117
 118         /// Sets all entries to x
 119         TProb<T> & fill(T x) {
 120             std::fill( _p.begin(), _p.end(), x );
 121             return *this;
 122         }
 123
 124         /// Draws all entries i.i.d. from a uniform distribution on [0,1)
 125         TProb<T> & randomize() {
 126             std::generate(_p.begin(), _p.end(), rnd_uniform);
 127             return *this;
 128         }
 129
 130         /// Returns length of the vector, i.e., the number of entries
 131         size_t size() const {
 132             return _p.size();
 133         }
 134
 135         /// Sets entries that are smaller than epsilon to 0
 136         TProb<T>& makeZero( Real epsilon ) {
 137             for( size_t i = 0; i < size(); i++ )
 138                 if( fabs(_p[i]) < epsilon )
 139                     _p[i] = 0;
 140             return *this;
 141         }
 142
 143         /// Set all entries to 1.0/size()
 144         TProb<T>& setUniform () {
 145             fill(1.0/size());
 146             return *this;
 147         }
 148
 149         /// Sets entries that are smaller than epsilon to epsilon
 150         TProb<T>& makePositive( Real epsilon ) {
 151             for( size_t i = 0; i < size(); i++ )
 152                 if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
 153                     _p[i] = epsilon;
 154             return *this;
 155         }
 156
 157         /// Multiplies each entry with scalar x
 158         TProb<T>& operator*= (T x) {
 159             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
 160             return *this;
 161         }
 162
 163         /// Returns product of *this with scalar x
 164         TProb<T> operator* (T x) const {
 165             TProb<T> prod( *this );
 166             prod *= x;
 167             return prod;
 168         }
 169
 170         /// Divides each entry by scalar x
 171         TProb<T>& operator/= (T x) {
 172             DAI_DEBASSERT( x != 0.0 );
 173             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
 174             return *this;
 175         }
 176
 177         /// Returns quotient of *this and scalar x
 178         TProb<T> operator/ (T x) const {
 179             TProb<T> quot( *this );
 180             quot /= x;
 181             return quot;
 182         }
 183
 184         /// Adds scalar x to each entry
 185         TProb<T>& operator+= (T x) {
 186             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
 187             return *this;
 188         }
 189
 190         /// Returns sum of *this and scalar x
 191         TProb<T> operator+ (T x) const {
 192             TProb<T> sum( *this );
 193             sum += x;
 194             return sum;
 195         }
 196
 197         /// Subtracts scalar x from each entry
 198         TProb<T>& operator-= (T x) {
 199             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
 200             return *this;
 201         }
 202
 203         /// Returns difference of *this and scalar x
 204         TProb<T> operator- (T x) const {
 205             TProb<T> diff( *this );
 206             diff -= x;
 207             return diff;
 208         }
 209
 210         /// Lexicographical comparison (sizes should be identical)
 211         bool operator<= (const TProb<T> & q) const {
 212             DAI_DEBASSERT( size() == q.size() );
 213             for( size_t i = 0; i < size(); i++ )
 214                 if( !(_p[i] <= q[i]) )
 215                     return false;
 216             return true;
 217         }
 218
 219         /// Pointwise multiplication with q (sizes should be identical)
 220         TProb<T>& operator*= (const TProb<T> & q) {
 221             DAI_DEBASSERT( size() == q.size() );
 222             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
 223             return *this;
 224         }
 225
 226         /// Return product of *this with q (sizes should be identical)
 227         TProb<T> operator* (const TProb<T> & q) const {
 228             DAI_DEBASSERT( size() == q.size() );
 229             TProb<T> prod( *this );
 230             prod *= q;
 231             return prod;
 232         }
 233
 234         /// Pointwise addition with q (sizes should be identical)
 235         TProb<T>& operator+= (const TProb<T> & q) {
 236             DAI_DEBASSERT( size() == q.size() );
 237             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
 238             return *this;
 239         }
 240
 241         /// Returns sum of *this and q (sizes should be identical)
 242         TProb<T> operator+ (const TProb<T> & q) const {
 243             DAI_DEBASSERT( size() == q.size() );
 244             TProb<T> sum( *this );
 245             sum += q;
 246             return sum;
 247         }
 248
 249         /// Pointwise subtraction of q (sizes should be identical)
 250         TProb<T>& operator-= (const TProb<T> & q) {
 251             DAI_DEBASSERT( size() == q.size() );
 252             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
 253             return *this;
 254         }
 255
 256         /// Return *this minus q (sizes should be identical)
 257         TProb<T> operator- (const TProb<T> & q) const {
 258             DAI_DEBASSERT( size() == q.size() );
 259             TProb<T> diff( *this );
 260             diff -= q;
 261             return diff;
 262         }
 263
 264         /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
 265         TProb<T>& operator/= (const TProb<T> & q) {
 266             DAI_DEBASSERT( size() == q.size() );
 267             for( size_t i = 0; i < size(); i++ ) {
 268                 if( q[i] == 0.0 )
 269                     _p[i] = 0.0;
 270                 else
 271                     _p[i] /= q[i];
 272             }
 273             return *this;
 274         }
 275
 276         /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
 277         TProb<T>& divide (const TProb<T> & q) {
 278             DAI_DEBASSERT( size() == q.size() );
 279             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
 280             return *this;
 281         }
 282
 283         /// Returns quotient of *this with q (sizes should be identical)
 284         TProb<T> operator/ (const TProb<T> & q) const {
 285             DAI_DEBASSERT( size() == q.size() );
 286             TProb<T> quot( *this );
 287             quot /= q;
 288             return quot;
 289         }
 290
 291         /// Returns pointwise inverse
 292         /** If zero==true; uses 1/0==0, otherwise 1/0==Inf.
 293          */
 294         TProb<T> inverse(bool zero=true) const {
 295             TProb<T> inv;
 296             inv._p.reserve( size() );
 297             if( zero )
 298                 for( size_t i = 0; i < size(); i++ )
 299                     inv._p.push_back( _p[i] == 0.0 ? 0.0 : 1.0 / _p[i] );
 300             else
 301                 for( size_t i = 0; i < size(); i++ )
 302                     inv._p.push_back( 1.0 / _p[i] );
 303             return inv;
 304         }
 305
 306         /// Raises entries to the power a
 307         TProb<T>& operator^= (Real a) {
 308             if( a != 1.0 )
 309                 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::ptr_fun<T, Real, T>(std::pow), a) );
 310             return *this;
 311         }
 312
 313         /// Returns *this raised to the power a
 314         TProb<T> operator^ (Real a) const {
 315             TProb<T> power(*this);
 316             power ^= a;
 317             return power;
 318         }
 319
 320         /// Returns pointwise signum
 321         TProb<T> sgn() const {
 322             TProb<T> x;
 323             x._p.reserve( size() );
 324             for( size_t i = 0; i < size(); i++ ) {
 325                 T s = 0;
 326                 if( _p[i] > 0 )
 327                     s = 1;
 328                 else if( _p[i] < 0 )
 329                     s = -1;
 330                 x._p.push_back( s );
 331             }
 332             return x;
 333         }
 334
 335         /// Returns pointwise absolute value
 336         TProb<T> abs() const {
 337             TProb<T> x;
 338             x._p.reserve( size() );
 339             for( size_t i = 0; i < size(); i++ )
 340                 x._p.push_back( _p[i] < 0 ? (-_p[i]) : _p[i] );
 341             return x;
 342         }
 343
 344         /// Applies exp pointwise
 345         const TProb<T>& takeExp() {
 346             std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::exp) );
 347             return *this;
 348         }
 349
 350         /// Applies log pointwise
 351         /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
 352          */
 353         const TProb<T>& takeLog(bool zero=false) {
 354             if( zero ) {
 355                 for( size_t i = 0; i < size(); i++ )
 356                    _p[i] = ( (_p[i] == 0.0) ? 0.0 : std::log( _p[i] ) );
 357             } else
 358                 std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::log) );
 359             return *this;
 360         }
 361
 362         /// Returns pointwise exp
 363         TProb<T> exp() const {
 364             TProb<T> e(*this);
 365             e.takeExp();
 366             return e;
 367         }
 368
 369         /// Returns pointwise log
 370         /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
 371          */
 372         TProb<T> log(bool zero=false) const {
 373             TProb<T> l(*this);
 374             l.takeLog(zero);
 375             return l;
 376         }
 377
 378         /// Returns sum of all entries
 379         T sum() const {
 380             T Z = std::accumulate( _p.begin(),  _p.end(), (T)0 );
 381             return Z;
 382         }
 383
 384         /// Return sum of absolute value of all entries
 385         T sumAbs() const {
 386             T s = 0;
 387             for( size_t i = 0; i < size(); i++ )
 388                 s += fabs( (Real) _p[i] );
 389             return s;
 390         }
 391
 392         /// Returns maximum absolute value of all entries
 393         T maxAbs() const {
 394             T Z = 0;
 395             for( size_t i = 0; i < size(); i++ ) {
 396                 Real mag = fabs( (Real) _p[i] );
 397                 if( mag > Z )
 398                     Z = mag;
 399             }
 400             return Z;
 401         }
 402
 403         /// Returns maximum value of all entries
 404         T max() const {
 405             T Z = *std::max_element( _p.begin(), _p.end() );
 406             return Z;
 407         }
 408
 409         /// Returns minimum value of all entries
 410         T min() const {
 411             T Z = *std::min_element( _p.begin(), _p.end() );
 412             return Z;
 413         }
 414
 415         /// Returns {arg,}maximum value
 416         std::pair<size_t,T> argmax() const {
 417             T max = _p[0];
 418             size_t arg = 0;
 419             for( size_t i = 1; i < size(); i++ ) {
 420               if( _p[i] > max ) {
 421                 max = _p[i];
 422                 arg = i;
 423               }
 424             }
 425             return std::make_pair(arg,max);
 426         }
 427
 428         /// Normalizes vector using the specified norm
 429         T normalize( NormType norm=NORMPROB ) {
 430             T Z = 0.0;
 431             if( norm == NORMPROB )
 432                 Z = sum();
 433             else if( norm == NORMLINF )
 434                 Z = maxAbs();
 435             if( Z == 0.0 )
 436                 DAI_THROW(NOT_NORMALIZABLE);
 437             else
 438                 *this /= Z;
 439             return Z;
 440         }
 441
 442         /// Returns normalized copy of *this, using the specified norm
 443         TProb<T> normalized( NormType norm = NORMPROB ) const {
 444             TProb<T> result(*this);
 445             result.normalize( norm );
 446             return result;
 447         }
 448
 449         /// Returns true if one or more entries are NaN
 450         bool hasNaNs() const {
 451             bool foundnan = false;
 452             for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
 453                 if( isnan( *x ) ) {
 454                     foundnan = true;
 455                     break;
 456                 }
 457             return foundnan;
 458         }
 459
 460         /// Returns true if one or more entries are negative
 461         bool hasNegatives() const {
 462             return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
 463         }
 464
 465         /// Returns entropy of *this
 466         Real entropy() const {
 467             Real S = 0.0;
 468             for( size_t i = 0; i < size(); i++ )
 469                 S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));
 470             return S;
 471         }
 472
 473         /// Returns a random index, according to the (normalized) distribution described by *this
 474         size_t draw() {
 475             double x = rnd_uniform() * sum();
 476             T s = 0;
 477             for( size_t i = 0; i < size(); i++ ) {
 478                 s += _p[i];
 479                 if( s > x )
 480                     return i;
 481             }
 482             return( size() - 1 );
 483         }
 484 };
 485
 486
 487 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
 488 /** \relates TProb
 489  */
 490 template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
 491     DAI_DEBASSERT( p.size() == q.size() );
 492     Real result = 0.0;
 493     switch( dt ) {
 494         case TProb<T>::DISTL1:
 495             for( size_t i = 0; i < p.size(); i++ )
 496                 result += fabs((Real)p[i] - (Real)q[i]);
 497             break;
 498
 499         case TProb<T>::DISTLINF:
 500             for( size_t i = 0; i < p.size(); i++ ) {
 501                 Real z = fabs((Real)p[i] - (Real)q[i]);
 502                 if( z > result )
 503                     result = z;
 504             }
 505             break;
 506
 507         case TProb<T>::DISTTV:
 508             for( size_t i = 0; i < p.size(); i++ )
 509                 result += fabs((Real)p[i] - (Real)q[i]);
 510             result *= 0.5;
 511             break;
 512
 513         case TProb<T>::DISTKL:
 514             for( size_t i = 0; i < p.size(); i++ ) {
 515                 if( p[i] != 0.0 )
 516                     result += p[i] * (std::log(p[i]) - std::log(q[i]));
 517             }
 518     }
 519     return result;
 520 }
 521
 522
 523 /// Writes a TProb<T> to an output stream
 524 /** \relates TProb
 525  */
 526 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
 527     os << "[";
 528     std::copy( P.p().begin(), P.p().end(), std::ostream_iterator<T>(os, " ") );
 529     os << "]";
 530     return os;
 531 }
 532
 533
 534 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
 535 /** \relates TProb
 536  */
 537 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
 538     assert( a.size() == b.size() );
 539     TProb<T> result( a.size() );
 540     for( size_t i = 0; i < a.size(); i++ )
 541         if( a[i] < b[i] )
 542             result[i] = a[i];
 543         else
 544             result[i] = b[i];
 545     return result;
 546 }
 547
 548
 549 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
 550 /** \relates TProb
 551  */
 552 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
 553     assert( a.size() == b.size() );
 554     TProb<T> result( a.size() );
 555     for( size_t i = 0; i < a.size(); i++ )
 556         if( a[i] > b[i] )
 557             result[i] = a[i];
 558         else
 559             result[i] = b[i];
 560     return result;
 561 }
 562
 563
 564 /// Represents a vector with entries of type Real.
 565 typedef TProb<Real> Prob;
 566
 567
 568 } // end of namespace dai
 569
 570
 571 #endif