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
  41
  42 namespace dai {
  43
  44
  45 /// Represents a vector with entries of type \a T.
  46 /** A TProb<T> is a std::vector<T> with an interface designed for dealing with probability mass functions.
  47  *  It is mainly used for representing measures on a finite outcome space, e.g., the probability
  48  *  distribution of a discrete random variable.
  49  *  \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
  50  */
  51 template <typename T> class TProb {
  52     private:
  53         /// The vector
  54         std::vector<T> _p;
  55
  56     public:
  57         /// Enumerates different ways of normalizing a probability measure.
  58         /**
  59          *  - NORMPROB means that the sum of all entries should be 1;
  60          *  - NORMLINF means that the maximum absolute value of all entries should be 1.
  61          */
  62         typedef enum { NORMPROB, NORMLINF } NormType;
  63         /// Enumerates different distance measures between probability measures.
  64         /**
  65          *  - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
  66          *  - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
  67          *  - DISTTV is the Total Variation distance;
  68          *  - DISTKL is the Kullback-Leibler distance.
  69          */
  70         typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
  71
  72         /// Default constructor
  73         TProb() : _p() {}
  74
  75         /// Construct uniform distribution over n outcomes, i.e., a vector of length n with each entry set to 1/n
  76         explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
  77
  78         /// Construct vector of length n with each entry set to p
  79         explicit TProb( size_t n, Real p ) : _p(n, (T)p) {}
  80
  81         /// Construct vector from a range
  82         /** \tparam Iterator Iterates over instances that can be cast to T.
  83          *  \param begin Points to first instance to be added.
  84          *  \param end Points just beyond last instance to be added.
  85          *  \param sizeHint For efficiency, the number of elements can be speficied by sizeHint.
  86          */
  87         template <typename Iterator>
  88         TProb( Iterator begin, Iterator end, size_t sizeHint=0 ) : _p() {
  89             _p.reserve( sizeHint );
  90             _p.insert( _p.begin(), begin, end );
  91         }
  92
  93         /// Returns a const reference to the vector
  94         const std::vector<T> & p() const { return _p; }
  95
  96         /// Returns a reference to the vector
  97         std::vector<T> & p() { return _p; }
  98
  99         /// Returns a copy of the i'th entry
 100         T operator[]( size_t i ) const {
 101 #ifdef DAI_DEBUG
 102             return _p.at(i);
 103 #else
 104             return _p[i];
 105 #endif
 106         }
 107
 108         /// Returns reference to the i'th entry
 109         T& operator[]( size_t i ) { return _p[i]; }
 110
 111         /// Sets all entries to x
 112         TProb<T> & fill(T x) {
 113             std::fill( _p.begin(), _p.end(), x );
 114             return *this;
 115         }
 116
 117         /// Draws all entries i.i.d. from a uniform distribution on [0,1)
 118         TProb<T> & randomize() {
 119             std::generate(_p.begin(), _p.end(), rnd_uniform);
 120             return *this;
 121         }
 122
 123         /// Returns length of the vector, i.e., the number of entries
 124         size_t size() const {
 125             return _p.size();
 126         }
 127
 128         /// Sets entries that are smaller than epsilon to 0
 129         TProb<T>& makeZero( Real epsilon ) {
 130             for( size_t i = 0; i < size(); i++ )
 131                 if( fabs(_p[i]) < epsilon )
 132                     _p[i] = 0;
 133             return *this;
 134         }
 135
 136         /// Sets entries that are smaller than epsilon to epsilon
 137         TProb<T>& makePositive( Real epsilon ) {
 138             for( size_t i = 0; i < size(); i++ )
 139                 if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
 140                     _p[i] = epsilon;
 141             return *this;
 142         }
 143
 144         /// Multiplies each entry with scalar x
 145         TProb<T>& operator*= (T x) {
 146             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
 147             return *this;
 148         }
 149
 150         /// Returns product of *this with scalar x
 151         TProb<T> operator* (T x) const {
 152             TProb<T> prod( *this );
 153             prod *= x;
 154             return prod;
 155         }
 156
 157         /// Divides each entry by scalar x
 158         TProb<T>& operator/= (T x) {
 159 #ifdef DAI_DEBUG
 160             assert( x != 0.0 );
 161 #endif
 162             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
 163             return *this;
 164         }
 165
 166         /// Returns quotient of *this and scalar x
 167         TProb<T> operator/ (T x) const {
 168             TProb<T> quot( *this );
 169             quot /= x;
 170             return quot;
 171         }
 172
 173         /// Adds scalar x to each entry
 174         TProb<T>& operator+= (T x) {
 175             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
 176             return *this;
 177         }
 178
 179         /// Returns sum of *this and scalar x
 180         TProb<T> operator+ (T x) const {
 181             TProb<T> sum( *this );
 182             sum += x;
 183             return sum;
 184         }
 185
 186         /// Subtracts scalar x from each entry
 187         TProb<T>& operator-= (T x) {
 188             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
 189             return *this;
 190         }
 191
 192         /// Returns difference of *this and scalar x
 193         TProb<T> operator- (T x) const {
 194             TProb<T> diff( *this );
 195             diff -= x;
 196             return diff;
 197         }
 198
 199         /// Lexicographical comparison (sizes should be identical)
 200         bool operator<= (const TProb<T> & q) const {
 201 #ifdef DAI_DEBUG
 202             assert( size() == q.size() );
 203 #endif
 204             for( size_t i = 0; i < size(); i++ )
 205                 if( !(_p[i] <= q[i]) )
 206                     return false;
 207             return true;
 208         }
 209
 210         /// Pointwise multiplication with q (sizes should be identical)
 211         TProb<T>& operator*= (const TProb<T> & q) {
 212 #ifdef DAI_DEBUG
 213             assert( size() == q.size() );
 214 #endif
 215             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
 216             return *this;
 217         }
 218
 219         /// Return product of *this with q (sizes should be identical)
 220         TProb<T> operator* (const TProb<T> & q) const {
 221 #ifdef DAI_DEBUG
 222             assert( size() == q.size() );
 223 #endif
 224             TProb<T> prod( *this );
 225             prod *= q;
 226             return prod;
 227         }
 228
 229         /// Pointwise addition with q (sizes should be identical)
 230         TProb<T>& operator+= (const TProb<T> & q) {
 231 #ifdef DAI_DEBUG
 232             assert( size() == q.size() );
 233 #endif
 234             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
 235             return *this;
 236         }
 237
 238         /// Returns sum of *this and q (sizes should be identical)
 239         TProb<T> operator+ (const TProb<T> & q) const {
 240 #ifdef DAI_DEBUG
 241             assert( size() == q.size() );
 242 #endif
 243             TProb<T> sum( *this );
 244             sum += q;
 245             return sum;
 246         }
 247
 248         /// Pointwise subtraction of q (sizes should be identical)
 249         TProb<T>& operator-= (const TProb<T> & q) {
 250 #ifdef DAI_DEBUG
 251             assert( size() == q.size() );
 252 #endif
 253             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
 254             return *this;
 255         }
 256
 257         /// Return *this minus q (sizes should be identical)
 258         TProb<T> operator- (const TProb<T> & q) const {
 259 #ifdef DAI_DEBUG
 260             assert( size() == q.size() );
 261 #endif
 262             TProb<T> diff( *this );
 263             diff -= q;
 264             return diff;
 265         }
 266
 267         /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
 268         TProb<T>& operator/= (const TProb<T> & q) {
 269 #ifdef DAI_DEBUG
 270             assert( size() == q.size() );
 271 #endif
 272             for( size_t i = 0; i < size(); i++ ) {
 273                 if( q[i] == 0.0 )
 274                     _p[i] = 0.0;
 275                 else
 276                     _p[i] /= q[i];
 277             }
 278             return *this;
 279         }
 280
 281         /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
 282         TProb<T>& divide (const TProb<T> & q) {
 283 #ifdef DAI_DEBUG
 284             assert( size() == q.size() );
 285 #endif
 286             std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
 287             return *this;
 288         }
 289
 290         /// Returns quotient of *this with q (sizes should be identical)
 291         TProb<T> operator/ (const TProb<T> & q) const {
 292 #ifdef DAI_DEBUG
 293             assert( size() == q.size() );
 294 #endif
 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 totalSum() const {
 389             T Z = std::accumulate( _p.begin(),  _p.end(), (T)0 );
 390             return Z;
 391         }
 392
 393         /// Returns maximum absolute value of all entries
 394         T maxAbs() const {
 395             T Z = 0;
 396             for( size_t i = 0; i < size(); i++ ) {
 397                 Real mag = fabs( (Real) _p[i] );
 398                 if( mag > Z )
 399                     Z = mag;
 400             }
 401             return Z;
 402         }
 403
 404         /// Returns maximum value of all entries
 405         T maxVal() const {
 406             T Z = *std::max_element( _p.begin(), _p.end() );
 407             return Z;
 408         }
 409
 410         /// Returns minimum value of all entries
 411         T minVal() const {
 412             T Z = *std::min_element( _p.begin(), _p.end() );
 413             return Z;
 414         }
 415
 416         /// Normalizes vector using the specified norm
 417         T normalize( NormType norm=NORMPROB ) {
 418             T Z = 0.0;
 419             if( norm == NORMPROB )
 420                 Z = totalSum();
 421             else if( norm == NORMLINF )
 422                 Z = maxAbs();
 423 #ifdef DAI_DEBUG
 424             assert( Z != 0.0 );
 425 #endif
 426             *this /= Z;
 427             return Z;
 428         }
 429
 430         /// Returns normalized copy of *this, using the specified norm
 431         TProb<T> normalized( NormType norm = NORMPROB ) const {
 432             TProb<T> result(*this);
 433             result.normalize( norm );
 434             return result;
 435         }
 436
 437         /// Returns true if one or more entries are NaN
 438         bool hasNaNs() const {
 439             return (std::find_if( _p.begin(), _p.end(), isnan ) != _p.end());
 440         }
 441
 442         /// Returns true if one or more entries are negative
 443         bool hasNegatives() const {
 444             return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
 445         }
 446
 447         /// Returns entropy of *this
 448         Real entropy() const {
 449             Real S = 0.0;
 450             for( size_t i = 0; i < size(); i++ )
 451                 S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));
 452             return S;
 453         }
 454
 455         /// Returns a random index, according to the (normalized) distribution described by *this
 456         size_t draw() {
 457             double x = rnd_uniform() * totalSum();
 458             T s = 0;
 459             for( size_t i = 0; i < size(); i++ ) {
 460                 s += _p[i];
 461                 if( s > x )
 462                     return i;
 463             }
 464             return( size() - 1 );
 465         }
 466 };
 467
 468
 469 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
 470 /** \relates TProb
 471  */
 472 template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
 473 #ifdef DAI_DEBUG
 474     assert( p.size() == q.size() );
 475 #endif
 476     Real result = 0.0;
 477     switch( dt ) {
 478         case TProb<T>::DISTL1:
 479             for( size_t i = 0; i < p.size(); i++ )
 480                 result += fabs((Real)p[i] - (Real)q[i]);
 481             break;
 482
 483         case TProb<T>::DISTLINF:
 484             for( size_t i = 0; i < p.size(); i++ ) {
 485                 Real z = fabs((Real)p[i] - (Real)q[i]);
 486                 if( z > result )
 487                     result = z;
 488             }
 489             break;
 490
 491         case TProb<T>::DISTTV:
 492             for( size_t i = 0; i < p.size(); i++ )
 493                 result += fabs((Real)p[i] - (Real)q[i]);
 494             result *= 0.5;
 495             break;
 496
 497         case TProb<T>::DISTKL:
 498             for( size_t i = 0; i < p.size(); i++ ) {
 499                 if( p[i] != 0.0 )
 500                     result += p[i] * (std::log(p[i]) - std::log(q[i]));
 501             }
 502     }
 503     return result;
 504 }
 505
 506
 507 /// Writes a TProb<T> to an output stream
 508 /** \relates TProb
 509  */
 510 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
 511     os << "[";
 512     std::copy( P.p().begin(), P.p().end(), std::ostream_iterator<T>(os, " ") );
 513     os << "]";
 514     return os;
 515 }
 516
 517
 518 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
 519 /** \relates TProb
 520  */
 521 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
 522     assert( a.size() == b.size() );
 523     TProb<T> result( a.size() );
 524     for( size_t i = 0; i < a.size(); i++ )
 525         if( a[i] < b[i] )
 526             result[i] = a[i];
 527         else
 528             result[i] = b[i];
 529     return result;
 530 }
 531
 532
 533 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
 534 /** \relates TProb
 535  */
 536 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
 537     assert( a.size() == b.size() );
 538     TProb<T> result( a.size() );
 539     for( size_t i = 0; i < a.size(); i++ )
 540         if( a[i] > b[i] )
 541             result[i] = a[i];
 542         else
 543             result[i] = b[i];
 544     return result;
 545 }
 546
 547
 548 /// Represents a vector with entries of type Real.
 549 typedef TProb<Real> Prob;
 550
 551
 552 } // end of namespace dai
 553
 554
 555 #endif