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