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