include/dai/prob.h

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